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1. Introduction. Bounds in the mean square sense for elastostatic boundary value 
problems have been found by Prager and Synge [1]f by a method in which they split the 
original boundary value problem into two relaxed problems such that the common solu- 
tion of the two relaxed problems is the solution of the original problem. They associate a 
set of six stress components with a point or vector in function space, and define the metric 
of the function space by a suitable choice of the scalar product. They show that the solu- 
tion of the original boundary value problem, interpreted as a vector in function space, 
must lie on a hypercircle. Since the metric of the function space is positive definite, the 
hypercircle is bounded, so that the location of the solution of the original boundary value 
problem on the hypercircle bounds it in the mean square sense. They further show that 
these bounds may be improved by introducing additional solutions of the relaxed prob- 
lems. 

In a later paper by Synge [2], the method of the hypercircle is applied to boundary 


value problems other than the elastostatic problems; these include the Dirichlet and 
Neumann problems. In that paper bounds are derived for the mean value of the solution 
in the neighborhood of a point. The method used differs from that of the present paper, 
but there are certain underlying ideas in common. 

In the present paper, we consider only the Dirichlet problem; we assume that the 
solution has already been located on a hypercircle in function space, using the general 
approach of the two preceding papers; we then proceed to determine bounds for the 


solution at a point. 

We consider the problem in Euclidean N-space Ey . Apart from any further signifi- 
cance of this generalization, it enables us to discuss simultaneously the two most interest- 
ing cases, N = 2, N = 3. 

In addition to determining bounds for the solution at a point, the same general method 
makes it possible to obtain bounds for the various derivatives of the solution. Bounds for 
the derivatives are considered separately in Sec. 4. 

The work up to this point is concerned specifically with the problem of determining 
bounds for the solution and its derivatives at a point interior to the domain of definition 
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of our problem. However, for a certain type of problem (e.g. the torsion problem), it is 
useful to have bounds for the derivatives at points on the boundary. In Sec. 5, bounds are 
obtained for the normal derivative at a point on a plane face. 

In Sec. 6, a numerical example is given to illustrate the practical application of the 
method. At a specified point, two bounds are obtained for the numerical value of the 
solution of the Dirichlet problem in a plane region. 

An alternative method, suggested by the work of Greenberg, is discussed in Sec. 7. 
This method is discussed for the two dimensional problem, but the reader may extend it 
to N-dimensions without difficulty. 

2. Notation and auxiliary theorem. Latin suffixes take the range 1, 2, --- , N, and 
the summation convention holds for repeated suffixes. The summation convention does 
not hold for Greek suffixes. The coordinates are rectangular cartesian, and differentiation 
with respect to a coordinate is indicated by a subscript following a comma (u,; = du/0z,). 

We shall denote by V an open domain in Euclidean N-space Ey , bounded by a closed 
surface B of (N — 1) dimensions: V may be simply or multiply connected. We shall have 
occasion to make use of an N-dimensional sphere with center at a general point P (con- 
tained in V) and with radius a. The N-dimensional interior of any such sphere will be 
denoted by v and its bounding (N — 1)-surface by b (Fig. 1). The unit normal 





Fig. 1. In Euclidean N-space, V is an arbitrary open domain bounded by an (N — 1)-dimensional surface 
B. The point P is surrounded by an N-dimensional sphere b of radius a and center at P: v is the volume 
enclosed by b. The unit normal n; is directed outward from V on B and outward from v on b. 


n,(t = 1,2, --- , N) will always be directed outward from V on B and outward from v on b. 
In all cases, integration will be denoted by a single sign and the range of integration 
indicated by the element of integration unless otherwise stated. 
We introduce the idea of function space, which we shall denote by F’, taking care to 
distinguish between a vector in F and a vector in the Euclidean N-space, Fy . The former 
will be denoted by heavy type, S; the latter will carry a subscript, u, . 
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We define a vector S in F to be any vector field p; in V + B. We require p, to satisfy 
certain differentiability conditions, namely, that V can be divided into a finite number 
of subregions such that p; is continuous and has continuous first derivatives in each such 
subregion. We further require that the normal component (p,;n,) of p; be continuous across 
the (NV — 1)-dimensional surface separating two such subregions. 

The scalar product of two vectors S and S’ will be denoted by S-S’ and defined by 


S-S’ = | ve dV, (2.1) 


where the volume element dV indicates that the integration extends throughout V. On 
putting S’ = S, we get 


isP?=S?=S8-S= | ve. dV. (2.2) 


This gives the metric in F and S = 0 implies p, = 0 throughout V. 
Synge [2] has established the fact that the extremity of the solution vector S lies on a 
hypercircle T in F and has expressed this fact by writing 
S=C+ Bj, (2.3) 
where the extremity of C is the center of T and RF its radius: C and R are considered to be 
known quantities for present purposes and are given explicitly by the following formulas: 


nt = 1 [se — Sosy? - Fisr-ny'] 


wel v=l1 


(2.4) 
1 mi ms 
Cas | s* - 2, ee) + 2 rset) |, 
w= v=1 
The vector S* is a known vector, and the vectors I), and If’ are two orthonormal sets of 
unit vectors which together form an orthonormal set of (m, + mz.) unit vectors, I, : 
I,-I, = 6,.,; (p,o0 = 1,2,---,m=m, + m,), (2.5) 


where 6,, is the Kronecker delta. The vector J is arbitrary except for the conditions ex- 
pressing that J is a unit vector, and that it is orthogonal to the known orthonormal set I, : 
J-J = |, J:1, = 0, {(p = 1, 2, “ites , m). (2.6) 

Let G be any vector in F and consider the scalar product S-G of G and the solution 
vector S. We determine maximum and minimum values for S-G as S ranges over the 


hypercircle I. 
We may represent the vector G as the sum of its projections on each of the unit 


vectors I, (9 = 1, 2, --- , m) and on the subspace of the hypercircle: that is 
G=MjJ,+ >N,I,, (2.7) 
p=l 


where 1/ J, (J, isa unit vector) is the vector projection of G on the plane of the hypercircle 
and 


N, = G-I,. (2.8) 
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To obtain M we square (2.7) and again use the orthogonality conditions; we obtain 


m 


M? = G’?— >> (G-I,)’. (2.9) 
p=l 
This value is always positive, so that J is real, and we shall understand by M the posi- 
tive square root of (2.9). Making use of (2.3), we may write 


SG-CG=RJ-G=R | 3, Pe 2 v,1, | =RMJ-J,. (2.10) 


Hence 


S-G —C-G| = RM | J-J.| < RM. (2.11) 


This may be summed up in the following theorem. 
Theorem I, Let G be an arbitrary vector in F and § the solution vector, the extremity of 
which is confined to lie on a hypercircle in F. Then bounds for the scalar product S-G are 


given by 
S-G —C-G| < RM, (2.12) 
Where R is the radius of the hypercircle and M is the positive square root of (2.9). 
3. Bounds for the solution at an interior point. The Dirichlet problem requires the 
determination of a function @ which is harmonic in a volume V and assumes assigned 
values on its surface B. These conditions may be written 


Ad = 0, (>)p = f, (3.1) 
where A is the Laplacian operator in /y and f is a function assigned on the (NV — 1)- 
dimensional bounding surface B; it will be assumed that f is piecewise continuous. The 
normal n, (7 |, 2, --- , N) to B will be assumed to be piecewise continuous and such 


that a solution ¢ exists which makes the various integrals occurring converge. 
It is well known that the fundamental solution of the Laplace equation is 


G = mr «., for N = 2, 
2—N (3.2) 
r » rf! 
= — flor: N > 3. 
(2 — N) = 
As a first step toward obtaining bounds for the solution ¢ at a point P with coordinates 
y; = ¢c; , we define in V a free Green’s vector, obtained by differentiation of the funda- 
mental solution: 
G, = 27°" forr > a, (4: 4 2. +++ , WV) 
(3.3) 
== forr <a. 
This formula applies for N = 2, 3, --- . The coordinates x; are relative to 
P(x; = y; — ¢c;), ris the distance measured from P(r* = x;x;), and a does not exceed the 
distance from P to the nearest point of B, so that the sphere r = a with center at P 


does not cut B. It should be noted that G,; defined by (3.3) is divergence-free, that is, 
G,,, = Nr™ — Naar"? = 0, (3.4) 


since x,x, = r°. This property is also evident from the fact that G; is the gradient of the 


fundamental harmonic function. 








of 
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We denote by G the vector in F corresponding to G; in V and associate the vector S 
in F with the vector field ¢,; in V, where ¢ is the solution of (3.1). We seek to bound @ 
at the point P. Noting that G; has no singularity in V, the scalar product S-G is evaluated 
by use of Green’s theorem: 


sG=[¢.,qv=|[ 9.G.aV 
_ (3.5) 
= [ oan, dB — | oan. db — | oG,;., dV, 
. Jy-s 
where db is an element of the sphere (r = a). The first surface integral is calculable 


since ¢ is known on B. The normal n; points outward from v on b and is equal to 2z;/r. 
Then by use of a mean value theorem for harmonic functions’, the integral over 6 is 
expressible in terms of the value of the function @ at the point P: 


| 6G n,; db = [ onc VY db = [ er’ v db 
(3.6) 


= a' z [ 6a _— Kyd@(P), 


where Ky is a constant depending only on the dimensionality of our space and is given 
explicitly by (3.9). 
We may now write (3.5) in the form: 


Ky(P) = | ¢G.n; dB — §-G. (3.7) 


This relation expresses the value of the solution at the point P in terms of a calculable 
quantity and the sealar product of the solution vector and the Green’s vector in the 
function space F, 

Theorem I of Section 2 provides bounds for the scalar product; so combining these 
results, we have bounds for the solution at the point P which are given by the following 
theorem. 

Theorem II. Let @ be a function which is harmonic in V and assumes assigned values on 
the boundary B of V. Let P be any interior point in V. Then bounds for ¢(P) are given by 
the inequality 

| sG.n, dB — C-G — Kyo(P)| < MR, (3.8) 
| « | 


! 


where n, is the unit normal outward from V on B. 
In (3.8), C, R are given by (2.4), M is the positive square root of (2.9), G; is given by 
(3.3), and 


; 2n*” . en 
Ky = TU?) : (K, = 2x, K; = 4n), 
f sail (>) 2 ’ (3.9) 


C-G = [ Cra” aV, 
“V-e 


10, D. Kellogg, Foundations of potential theory, Berlin, 1929, pp. 223-224. 
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where C’;; denotes the vector field in V corresponding to the vector C in F whose extremity 
lies at the center of the hypercircle. 

4. Bounds for the derivatives of the solution at an interior point. ‘To obtain bounds 
for the first derivatives, we define a new Green’s vector G;”’ which is derived from the 
Green’s vector of Section 3 by differentiation of equation (3.3) with respect to x, : 


y(z a \ 4 \ 4 ¥ 
Gr =t,9" —Nia27 , tory > a, 


(4.1) 


0, forr < a, (@,p = 1,2,---, WN). 


It is of course obvious that this vector is divergence-free since it is the second deriva- 
tive of the fundamental harmonic function. 

We denote by G’”’ the vector in the function space F corresponding to G;”’ in V and 
evaluate the scalar product S-G’”, making use of Green’s theorem: 


S-G = | 6 G; dV = | Gy dV 
ny oy 
(4.2) 


= [ ¢G{"'n, dB — [ oG?'n, db — | 


. /y- 


eG”, aV. 
: 

The integral over B is calculable since ¢ is known on B, and the last integral vanishes 
since G;”’ is divergence-free. We evaluate the remaining integral, making use of Green’s 
theorem and another mean value theorem for harmonic functions’. We observe that 


n, = 2;/r and N;X, = NX; ON b. Then 


| 6G.” n; db = a“ on, db — Na“ it ox ,x,n; db 


ny 


= (1 — N)a” | ¢,,dv = (1 — N)Ly @¢,,(P), 


where Ly depends only on the dimensionality of the space and is given explicitly by 

s : vr 

Ly = Ky : Ses wae (4.4) 
. r(N/2 + 1) 

Equation (4.2) is now in the form: 


(1 — N)Ly¢.,(P) = | oG;"n; dB — §-G’”’. (4.5) 

This relation gives an expression for the value of the first derivative of the solution at the 
point P in terms of a calculable integral and the scalar product of the solution vector and a 
Green’s vector in the function space F. 

The result (4.5) when combined with Theorem I gives bounds for the first derivatives 
of the solution ¢ at the point P. Thus, we have established the following theorem. 

Theorem III. Let be a function which ts harmonic in V and assumes assigned values on 
the boundary B of V. Then bounds for the first derivatives of the solution at a point P con- 
tained in V are given by the inequality 


*Loc. cit. pp. 223-224. 
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° | 
| {G2??n, dB — C-G™ + (N — I)Ly¢,,(P) | < MR, (4.6) 
d | 

where n,; represents the unit normal outward from V on B. 

The value of Ly is given by (4.4) and M"” is obtained from (2.9) by replacing G by 
. 

A process similar to the preceding gives bounds for the higher order derivatives of ¢ 
at a point P interior to V. The appropriate vector field in V is obtained in each case by 
differentiation of (3.3). The variables with respect to which (3.3) is to be differentiated are 
determined by the particular derivative of ¢ for which bounds are being sought. For exam- 
ple, if bounds are required for ¢,,, , then the vector field is obtained for this case by 
differentiation of (3.3) with respect to x, and x, . The resulting vector field is then made 
to correspond to a Green’s vector in the function space F and the scalar product of this 
Green’s vector and the solution vector is computed. Combining the results of this compu- 
tation with Theorem I gives bounds for the derivative under consideration. Such results 
may be stated in a manner analogous to Theorem III with appropriate changes. Thus, if 
bounds were required for ¢,,, the inequality (4.6) must be replaced by 

N? = (va) , 
VN+2 Ly@,»¢P) | cn (4.7) 
with the obvious meaning given to the changed notation. 


5. Bounds for the normal derivative at a point on a plane boundary. The preceding 
work is concerned with the problem of determining bounds for the solution and its 


| sG2n, dB — C-6° — 


> 











Fic. 2. Ad = 0 in V, (¢)g = f, a given function. The origin of coordinates is at P and the 2x;-axis is per- 
pendicular to the plane portion of the bounding surface B. The hemispherical domain surrounding P is 
bounded by 6 + c, where c = B — B’. 


derivatives at a point interior to the domain V. However in certain types of problems, 
it is useful to have bounds for the normal derivative at a point on the boundary. In this 
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section, we confine our attention to a domain V bounded by a closed surface B such that 
at least one portion of B is a plane (Fig. 2). We shall establish bounds for the normal 
derivative of the solution at a point P on this plane portion of the bounding surface. 

A system of rectangular cartesian coordinates is introduced such that the origin is at P 
and the positive x,-axis is perpendicular to B and is directed outward from V. Thus the 
derivative ¢,, is the normal derivative of ¢ at P. The point P is then surrounded by a 
hemispherical domain v of radius a lying entirely in V. The surface of the hemisphere is 
denoted by b and that portion of B which caps b is called c. We use B’ to denote B — c. 

We define a Green’s vector G‘”’ in the function space F corresponding to the vector 
field in V obtained from (4.1) by setting p = 1, and compute the scalar product of G™ 
and the solution vector S. However after applying Green’s theorem, the surface integrals 
over B and the boundary of v have the range of integration c in common so that the scalar 
product takes the form 


. . 


S-G = | 9G(0n, dB’ — | 6G{n, db. (5.1) 
As in the previous instances, the first integral is calculable since ¢ = f on B’, but we no 
longer have the mean value property of harmonic functions to aid in the evaluation of the 
last integral since it is an integral over the hemispherical surface b. However, this last 
integral can be evaluated by other means to obtain the result 


[ oG ». db=a . [2 ' dz [ F(x) dx + (N — Ii)a 5 [ 7 de 


where 


F(x) = (N — 1)2' | (b.22 + 6,53 + +++ + o.ny) de. (5.3) 


Since ¢ = f onc and the 2,-axis is perpendicular to c, the integrand of (5.3) is a known 
function, hence F is a known function. We note that this evaluation of the last integral 
of (5.1) requires that the boundary function f and its first and second derivatives be 
piecewise continuous and such that their integrals exist. 

Upon substitution from (5.2) into (5.1), an explicit expression for S-G'"’ is obtained 
in terms of calculable integrals and the normal derivative of ¢ at the point P. Then com- 
bining this result with Theorem I, we obtain the following theorem. 

Theorem IV. Leto be a function which is harmonic in V and assumes assigned values on 
the boundary B of V. It is further required that B contains a plane portion. Then bounds for 
the normal derivative of ¢ at a point P on this plane portion are given by the inequality 


- ra “ 


| {G0?n, dB’ — C-G" — a” | 2 dz | F(x) dx 
~ B’ “0 “0 
(5.4) 
y — F ] , 0¢(P) | ()p 
—(N - 1a } fde—5(1 — N)Ly “S| < MR, 


where n,; is the unit outward normal. 
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It is interesting to note that Theorem IV gives an exact value in the following case. 
Let ¢ be harmonic in the domain V, where V is a volume in Euclidean N-space bounded 
by a hemispherical surface B’ of radius a and a diametral plane c. The point P is taken at 
the center of curvature of B’ (Fig. 3). We assign ¢ the value zero on ¢ and unity on B’ 


Bt 








Fia. 3. Ad = Oin V, ()g- = 1, )- = 0. The normal derivative of ¢ at the point P is computed for the 
special case in which V consists of a domain bounded by a hemispherical surface and a diametral plane in 
Euclidean N-space. 


and seek the value of the normal derivative of ¢ at the point P. We may take the volume v 
to coincide with V, so that obviously M“™ and C-G“” of (5.4) are zero. Further, since 
@ = 0 one, it follows that F(a) = 0 and that the integral of f over c is also zero. The 
remaining integral over B’, when evaluated, gives 





sii (1 — N)r-?? ” 
G}' = . oO 
I, JG, n, dB’ = “orw — 1/2] (5.5) 
Substitution of (5.5) into (5.4) gives 
(N — 1)Ly- - 
¢@,(P) = a . (5.6) 
For N = 2, N = 3, this gives respectively 
¢,(P) = 2/(az), (N = 2), 
(5.7) 
¢(P) = 3/(2a), (N =3). 
This quantity represents the intensity at the center of a circular (NV = 2) or spherical 
(N = 3) condenser which is split symmetrically and the two halves raised to potentials 


+1 and —1. 

6. A numerical example. As an illustration of the use of the method contained in the 
preceding sections, we consider the following problem: let the open domain V be the area 
in the Euclidean plane (V = 2) bounded by B, x = +1, 2” + y’ = 2 (Fig. 4). We seek 
bounds for the function y which satisfies the conditions 


Ay=0, We =f=3@'+y’). (6.1) 


In particular, we establish numerical bounds for the function y at the origin (the point P 
in the notation of the preceding work). 
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i] 
nN 
i) 


We use the explicit (x, y) notation and the implicit notation interchangeably, depend- 


ing on which is the more convenient. 
We first find a vector field p* and a scalar field y* in V satisfying the conditions 


pPr= vi, (2 =f. (6.2) 
Secondly, we find a vector field p/’ in V satisfying 
p,’; =0 (6.3) 
and no boundary conditions. The common solution of (6.2) and (6.3) is the solution of 
(6.1). We can obtain solutions of (6.3) by means of skew differentiation: p/’ = €;;x’/ . 
y 


AR 
y 


Fig. 4. AY = Oin V, (¥)g = 3(x? + y?). The domain V is located in the Euclidean plane and is bounded 
by B: (x? + y? — 2)(z? — 1) = 0. The point P is taken at the origin and v consists of the interior of the 














(7 





unit circle. 


For the purpose of improving the bounds, we use a set of vectors of the class (6.3) 
and introduce another class of vector fields p/ in V satisfying the conditions 


p=wWv; (d’), =, (6.4) 


Thus our problem is that of selection of the three types of functions from which we 
can derive the vector fields p* , p;’, p: satisfying the appropriate conditions. The func- 
tions selected for use in our particular example are given in Table I. This table shows the 
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vector in F’, the vector field in V to which it corresponds, and the particular function from 
which the vector field is derived. The functions given in Table I are used to obtain only 
two approximations for the solution function y. It is, of course, possible to improve the 
approximation last obtained by the selection of additional functions. 











i _ TABLE a 
Vector , | Function from which 
in F Vector field in V vector field is derived 
S* | pt =a | 
y* = (2° + y’)/2 
pi=y | 





Vay = (2 — I(e’? + y’ — 2) 





S! Ploy, = 2u(Bx* + Qx*y? — 6x? — y? + 2) | 








| Vio) = a°(x — 1)(2” + y’ — 2) 
Dloyo = 2x°y(x° — 1) 
S Dis, = Qry*(2Qx* + y* — 3) 
| Vian = y*(a® — Ia? + y? — 2) 
Piso = Qy(x” — 1)(x” + 2y’ — 2) | 
Ss’ lov, =2 
ta = ty 





S2’ pin =x — Say’ 
” 3 3 
| xt = zy — ay 

| 


Dine = y° — 3x°y 








To determine what computations are required, we examine the quantities contained 
in the inequality (3.8) of Theorem II, Section 3; namely, for N = 2, 


| [ {Gn;,dB — C-G — 2ry¥(P) | < MR. (6.5) 


It is convenient to have all the scalar products available before computing the various 
quantities present in this inequality. These scalar products have been tabulated in a 
triangular array in Table II. The number in Table II which represents the scalar product 
of any two vectors S, S’ is located at the intersection of row S and column S’ (or vice 


versa). 
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TABLE IT. Scalar products. 
s+ | 4.474926 
s; |—s.633033 | 21.474728 
s; |—1412487 | 1.458289 | 3.253382 | 
s; |—3.170796 | 6.953487 | 395429 | 11.049406 | 
sy |-1.333333/ 0 |o | 0 | 4474926 
Si’ |— .533333 0 : 0 : 0 — 2.133333 | 7.654613 
G | 2.000000 |—3.141593 |—.785398 |—2.356194 |— .858407 |—.050148 |1.493808 | | 
ieee SS es ee, eee ees = a Sa | 
s* Si S; S; Sf’ S2’ 


- | 
We next orthonormalize the vectors S/ , S{’. It is inherent in the method that 
S/-S/’ = 0 for all 7, 7 (this may easily be checked by integration), so it is only necessary 
to require that 
S/. S/ 
(6.6) 
S7? F ag — 6 


A systematic method of orthonormalization has been given by Peach [7]. We apply his 
A syst t thod of ortl 1alization has beer en by Peach [7]. W ly | 
method to the two preceding classes of vectors S/ , S/’ to obtain the two orthonormal 


TABLE III. Scalar products. 














I; I; I; Ii’ I)’ 
gx | —1 863019 | 465193 | —.133665 | —.630298 | 453704 
~~~ | 677944 | s22t03 | —.456269 | —405789 | 178320 | 

(S*-I/)? = 3.470839 (G-Ij) = .459608 

(S*-I3)? = .216405 (G-Ii)” = .103750 

(S*-1)? = .017866 (G-I)* = .208181 

(S*-1/’)? = 397275 (G-Ii’)? = .164665 

(S*-14’)? = .205929 (G- Ii’)? = .031801 
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classes I , I/’ (obviously, I{-I/’ = 0, for all 7, 7). When these computations have been 
carried out, we find: 


Ii = .215792 Si, 
I; = —.038235 S{ + .563047 S; , 
I; = —.109735 S} + .008205 S: + .337178 S;, 
(6.7) 
Ij’ = 472723 Si’ , 
Ii’ = __-.185041 Si’ + .388146 Si’ . 


We also need the scalar products of S* with each I’, I’’, and the square of all such 
products. These are easily obtained from (6.7) and Table II: the results are given in 
Table III. 

The approximations for R® are given by direct substitution from Table III into the 
first of (2.4). The first two approximations for R’ are as follows: 


9 


R, = .151703, (approximation based on S*, Ij , Ij’), 
(6.8) 


Il 


R; = .041653, (approximation based on all vectors given above). 

It should be pointed out that the above quantities are independent of the position of the 
point at which we are seeking to bound y. In fact, both the radius of the hypercircle and 
the position of its center depend only on the domain V and on the functions from which 
the vector fields are obtained. Consequently, we would write down the first two approxi- 
mations for C if we so desired. However, we are not directly interested in C, so we proceed 
to the computation of C-G and M. 

We use the Green’s vector defined by (3.3) for N = 2 and a = 1; that is, 


© ¢ 


G, = zr 
forr > 1, 
G. = yr” (6.9) 


G, = G, = 0, forr < 1. 


The necessary scalar products of G with each of the other vectors have been com- 
puted and are included in the appropriate tables. We then form the scalar product C-G 
by dotting G into the second of (2.4). Then direct substitution from Tables II and III 
gives the first two numerical approximations for C-G. Substitution from Table III into 
(2.9) gives the approximations for M*. These results are as follows: 


My = .869535, M, 


.932488, (first approximation), 


I 


G-C, = .496371, (first approximation), 
(6.10) 


M? = .525803, M, = .725122, (second approximation), 


G-C, = .431422, (second approximation). 
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The line integral over the boundary B of the domain V is found to be 


+ 


| fGn;dB =2+-. (6.11) 
We may now substitute from (6.11), (6.10), (6.8) into (6.5) to obtain bounds for the 
function y at the point P(0, 0). This gives 


682 < ¥(P) < .797, 


(6.12) 


(26 S WP) < .773, 


as the first and second approximations respectively. 

7. An alternative method. The reading of the work of H. J. Greenberg [8] on the 
same subject suggested an alternative method of attack on the problem of obtaining 
bounds at an interior point for the solution of the Dirichlet problem. The method is out- 
lined below for the special case N = 2, but may obviously be extended to the N-dimen- 
sional case. 

We let ¢ represent the solution we seek to bound at the point P and S the solution 
vector in function space corresponding to ¢. We use two functions: the first is the free 
Green’s function G, and the second any function H which is regular (the function and its 
first derivatives continuous) throughout V and is equal to G on the boundary B of V. 
Thus we have: 


G = logr, 


(H), = (G)z . 


We consider the scalar product S-G of the solution vector S and the Green’s vector G. 
We cannot apply the divergence theorem to S-G directly so we surround P with a small 
sphere ) with center at P and radius a and remove the interior v of b. Then applying the 
divergence theorem to S-G for the region V — v in both possible ways and making use 


of the fact that both ¢ and G are harmonic, we have 


i 
Pe | P “RB 


[ ¢.G,,dV = [ ¢G..n; dB — [ oG..n; db 


a; 


= [ Gon; dB — | Go, in, db. 


© F 


Upon taking the limit as the radius a goes to zero, we obtain the following well known 
expression for the value of ¢ at point P: 


2ng(P) = | $G,.n, dB — | Go,.n, dB. 


vf vE 


~J 


The first integral is calculable since ¢ is known on B. The essential point of the method 
concerns the second integral. Since the regular function H is equal to G on B, we replace 
G by H and again apply the divergence theorem. Since ¢ is harmonic, this gives, 


| Go.n. dB = | Ho, dB = | H.,.dV = ELS. (7.4) 
B “B “V 


. 








n 
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Thus equation (7.3) may be written: 
ong(P) = | $G,n, dB — HS. (7.5) 
“B 


Now if the extremity of S is confined to a known hypercircle (cf. Section 2), we know that 
H-S is bounded above and below. Hence (7.5) gives bounds for ¢(P). 

It will be recalled that in the work of Section 2 the norm of the Green’s vector | G | 
enters into the bounds for G-S via the number M (see equation (2.9)); and in order to 
insure the existence of | G |, it was necessary to define our vector field to be zero in a region 
v containing P. This difficulty does not now arise, because it is | H |, not | G |, that enters, 
and | H | is finite since H is regular. 

The question of practical importance now is whether we can find the function H 
which has the required properties. This question has already been answered by Greenberg 
since one of his functions has these precise properties. The function referred to is defined 
in the following way: 


H = logr, forr > a, 
(7.6) 
=ear+art+e, forr <a, 


where the coefficients ¢, , c2 , ¢; are determined so that H has continuous second deriva- 


tives across b, that is: 


Cc; => —a™*/2, 
C, = 2a, (7.7) 
c, = loga — 3/2. 


But this is not the only way in which H may be defined, at least for convex regions V. 
For any point P interior to V, the distance r from P to any point on the boundary B may 


be defined by a single valued function of a polar angle @ at P; explicitly, 


r = p(8), 0 < @< 2r. (7.8) 
Then, if we define H by the relation 
r’ logr 
7 = ——; + h, 7.9 
(e(8))” ” ss 


where h is harmonic in V and vanishes on B, we have a function such that (H), = (@)z 
and has the required regularity at P. 

In conclusion, the author acknowledges his indebtedness to Professor John L. Synge 
under whose supervision this paper was prepared. 
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THE ELECTROMAGNETIC FIELD PRODUCED BY A HELIX* 


BY 
R. S. PHILLIPS 
University of Southern California 


1. Introduction. In this paper’ we are concerned with finding the electromagnetic 
field produced by an infinitely-thin perfectly-conducting helix subject to a mono- 
chromatic excitation. Because of the use of the helix in an ultra-high frequency amplifier 
known as the traveling-wave tube [2], there has been in recent years considerable interest 
in this problem. Earlier papers dealing with the field due to a helix have been based on 
a model consisting of a circular cylinder which conducts only in a given helical direction.” 
[1], [2], [3]. The symmetry of this idealization corresponds to that of cylindrical co- 
ordinates so that such a model yields to the usual separation of variables technique. 
Our approach is quite different. 

We start off by assuming that the effect of a monochromatic electromagnetic field 
would be to induce a sinusoidal electric current on the helical conductor; and that this 
current travels along the helix with a real propagation constant @ in the axial direction.* 
The electromagnetic field accompanying such a current may then be represented by 
means of the retarded vector potential IT as [5, p. 165] 


M(x, y,2,) = A | (l/r) ds, (1.1) 


where [7] is the current vector at a point (&, n, ¢) of the helix at time (¢ — r/c), c being 
the free space wave velocity, r = [(x — &)” + (y — »)? + (2 — §)’]'”, s the are length 
along the helix, and A a scalar factor. 

It follows that 


E=V(V-M+k0, H = —iw VXI, (1.2) 


where k° = euw , w being the radial frequency of the field. We then attempt to determine 
8 so as to make the electric field along the prefectly conducting helix “essentially”’ zero, 
and thus validate our initial assumption. It will turn out that the 8 so determined corre- 
sponds to a phase velocity along the helix itself equal to c. 

This then is our program. The principal difficulty lies in the singular nature of the 
potential function II. In Sec. 2 we show that the integral defining II converges; and in 
addition we obtain an evaluation for II in terms of Bessel functions. In fact we establish 


€ as mae R) dg = > pe Ar e"*1,(T,.a)K,(T,p), (p > a), (1.3) 


*Received June 20, 1949. 

1The work on this paper was done under funds granted to New York University by the Air Force 
Cambridge Research Laboratories, under contract W28-097-ac-170. 

2By helical direction we mean the direction of a vector field which in cylindrical coordinates has the 
components v, = 0, v, = a,v, = a. 

3That is, the amplitude varies as cos (8z — wt) where z is measured in the axial direction and w is 
the radial frequency. 
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where 
2) 2 9 ee, 1/ 
+a +p — 2ap cos (y — £/a)] 


and 

Tr, = [(6 — n/a)? — k’]'”. 
Here, J, and K, are modified Bessel and Hankel functions of nth order. If for a partic- 
ular n, I, is imaginary, then the corresponding term in (1.3) is to be replaced by 
Qn ie"? J,(y,a)H\” (y,p), where y, = [k° — (8 — n/a)’]'”. 

In Sec. 3 we prove that II may be differentiated either by differentiating under the 
integral sign or by differentiating its series expansion termwise. We then show that II 
satisfies the homogeneous wave equation at all points not on the helix. Also in Sec. 3 
explicit expressions for the electromagnetic field are obtained. Finally in Sec. 4 the 
propagation constant 8 is determined so that the boundary conditions for the helix are 
satisfied. Since the electric field is singular on the wire, 8 is determined so that the ratio 
of the electric field in the helical direction for the proper 8 to that of the electric field 
in the helical direction for any other 8 approaches zero as the helix is approached radially. 

In cylindrical coordinates, the equation for the helix is 


p = a4, z= ag. (1.4) 
The assumed current along the helix is the real part of 


i, =0 


] 


i, = Tala? + a’)"” exp [i(Bz — w)], (1.5) 


¢ 


i, = Tala’ + a’ |” exp [7(Bz — wf)]. 


The retarded potential at the point (p, ¢, z) and the time ¢ is the real part of 


p 


ll, = aF | r-* exp (i8¢ + zkr) sin (g — ¢£/a) df, 
ll, = aF | r* exp (i8¢ + tkr) cos (¢ — ¢/a) dé, (1.6) 


Il. =aF | r* exp (i8¢ + tkr) dé, 


where 

r= [2 — §)? +a’ + p” — 2ap cos (y — £/a)]'” 
and 

F = —Ie~**'/A4riwex. 


If we replace sines and cosines by their exponential representations, we see that all 
of the integrals in (1.6) are essential of the type found in II, . The notation of the next 


- 


section will be somewhat simplified if we make a change of variables ¢’ = ¢ — z in these 
integrals. Then 
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[ee /r) dg = | efF'(e' fr’) dt", (1.7) 


where 
= [t? + a + p* — 2apcos[(¢y — z/a) — ¢’/al]'”. 


2. Evaluation of the integral. In this section we shall evaluate the integral 


I= | e*¥(o**® /R) de, (2.1) 


where 
R = [¢ + a -e p = 2ap cos (yg ae ¢/a))'”. 


We restrict our consideration to real 8 since the integral will clearly not converge for 
complex £8. 

The usual method of handling improper integrals, namely to transform the path of 
integration, is not suitable here because R has an infinite set of branch points. These 
are the solutions of the transcendental equation R’ = 0. The device which we shall use 
in evaluating (2.1) consists of replacing the integral along the helix by a double integral 
over the enveloping cylinder of the helix; this leads to a Bessel function solution. 

We define 


0 L Qe 
uL,r =| as | doe*(e**/R)P(r, ¢ — ¢/a), (2.3) 
v-L “0 
where P(r, c) is the Poisson kernel [7, p. 51] 


P(r, 0) = 2 a oF 


Q2r1+r — 2rcose oF coe 





and 


R = [¢? + a” + p’ — 2ap cos (y — o))"” 


For 0 < r < 1, L <@, p > a, the double integral is absolutely integrable and hence 
the order of integration may be interchanged. 
By well known properties of the Poisson kernel 


lim | da e'#(e*" /R)P(r, o — t/a) = e“*F(e'** /R) 


r-1— “0 
uniformly for —L < ¢ < L. Hence 


I = lim lim J(L, r) (2.4) 
Rew 


©» r-l— 


if this limit exists. 
We wish to show not only that the limit of (2.4) exists but that the order of this 
double limit may be interchanged. In order to accomplish this we shall prove 


a) lim J(Z, r) exists for each r (0 < r < 1) 
Lo 


b) lim J(L, r) exists uniformly in L for (0 < L <@), 
r +1 








232 R. S. PHILLIPS [Vol. VIII, No. 3 


Since the Fourier expansion of P(r, «) converges uniformly in o for each r < 1, it 
follows that we can interchange the order of integration and summation and so express 
I(L, r) as 


x ol o25 
I(L,r) = > ri | dé | do R~* exp (7B + tkM) exp [in(o — ¢/a)]. (2.5) 
r s) “ L “0 
In order to establish (a) we first examine the individual terms of this sum. For the 


case (8 — n/a) > k, we make the substitution 


c [a” + p’ — 2ap cos (gy — o)|"’ sinh (y — 4), 


II 


II 


HR [a> + p — 2ap cos (gy — a)]"”* cosh (y — 5), 
where 
cosh 6 = (8 — n/a)/[(8 — n/a)* — k*}"”, 
sinh 6 = k/[(8 — n/a)’ — k’]'”. 
Then 
(8B — n/a)i + kM = SG sinh y, 


where we have set 


© = [a + p — 2ap cos (yg — a)]' (8 — n la)? — Ry”. (2.6) 
Hence, 

2 L als 

| dt (e**/®) exp [7(8 — n/a)t] = | dy exp (¢@ sinh y), 

vag J =t 


where /, and /, are defined by 
L = [a® + p* — 2ap cos (y — o)]'”’ sinh (1, + 4), 


e’ 


L = [a? + p* — 2ap cos (y — a)]'” sinh (I, — 4), 
5 being defined as above. We now transform the path of this integral to the line y = i/2. 
Then because of the fact that the integral along the line joining (/, 0) and (1, 7r/2) goes 
to zero as! — +, we get 


| dt (e**/®) exp [i(8 — n/a)f] = | dy exp (—® cosh yp) = 2K,(®). (2.7) 


The latter equality can be found in [6, p. 182, (7)]. 
It is clear from the above argument that 


ol 


2K,(®) — | . dt (e**/R) exp [i(8 — n/a)t] = | dy exp (t®@sinh y), 
J—L oF? 


where P is the sum of the two paths indicated on Fig. 1. 











bo 
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~ 


For p > a, the integral 


| dy exp (t® sinh y) 
Jp 


converges to zero uniformly in o as l, , 1, +. This incidentally shows that the order 
of integration in 


| do | dt e'*e**/R) exp [in(o — £/a) 


may be interchanged. We can, however, easily obtain an upper bound in the value of 
this integral over P. In fact the integral from (J, 0) to (1, ri/2) is bounded by 2/(# cosh 1); 














P p , 
re " incre mi/2 
tb, l, 

Fig. 1 


whereas the integral from (I, ri/2) to (@, ri/2) is bounded by e” *'/®. Now for sufficiently 
large | n |, fixed 8, and fixed p > a, ® > c|n|. Hence in this case 


.2" eL 


| da K,(®) e'”” — | d¢ | da e**(e**/) exp [in(o — £/a)] 


| 
te | Pee * | 


< = (2.8) 
| n | 


A similar argument shows that (2.7) and (2.8) hold equally well for (@ — n/2) < —k. 
For the case (8 — n/2)? < k’, let 


¢ = [a’ + p* — 2ap cos (y — o)]'” sinh (y — 8), 


R [a> + p’ — 2ap cos (v — a)]'”” cosh (y — 8), 


where now 
sinh 6’ = (8 — n/a)/[k? — (8 — n/a)*}”, 
cosh 6’ = k/[k? — (B — n/a)*}'”. 
For 
®’ = [a? + p® — 2ap cos (y — o)]'"[k’? — (8 — n/a)*]'” 
we have, then, 


(8 — n/a)i + kR = & cosh y 
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so that 


al pi’s 
| dé (e*"/R) exp [7(8 — n/a)gé] = | dy exp (7®’ cosh y). 
J-] J’, 

The difference between this integral and the integral from (— © — ix/2) to (—7m/2) 


to (ir/2) to (~ + in/2) is the integral over the path P’ shown in Fig. 2. 














p’ 
——ee i 77/2 
-1; 4 
: ly 
17/2 
id 
Fig. 2 


For p > a, the integral [p. dy exp (i&’ cosh y) converges to zero uniformly in a as |, , 


l, ~o. Hence 
| le (e"*"/R) exp [7(8 — n/a)ft] = | dy exp (i®’ cosh y) 
Aa a (2.9) 
= riH,'’(%’) 


(see [6, p. 180, (8)]) and 


| dt | da e'(e*" /R exp [in(o — ¢/a)] = | do H‘(®’)e'"’. (2.10) 


“0 


Now by the addition theorem for Bessel functions [6, p. 361] for p > a 


H;(a[a’ + p* — 2ap cos 6]'”*) = Zz, J,,(va)H (xp)e*”’, 


K,(zla° + p — 2ap cos 6]’") = Zz T,,(va) K,,(xp)e"”’, 


and these series converge uniformly in 6. Substituting these series in the above integrals 


we obtain for p > a 
129 


bout sBty ak : 
= | dé | do ¢ at St) exp [in(o cand a) | 


ai « x “9 


21/2 


27,([(8 — n/a)” — k*]'’a)K,(((8 — n/a)” — k*]'’*p) 


(2.11) 


= mid,((k? — (8 — n/a)*|'a)Hi" ({k’ — (8 — n/a)*]'”’p) 


if (8B — n/a)’ < k’. 
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Now it can be shown [4, corollary 2.2] that J,(x)K,(x) is a monotonic decreasing 
function of x which assumes the value 1/| 2n| at x = 0. Since, in addition, K,(x) is 
monotonic decreasing in 2, it follows that 


I,(T,a)K.(Tp) < 1/| 2n | 


for p > a. Combining this with the above result and the inequality (2.8), we obtain 
for fixed 8 and fixed p > a 


| L Qr 
/ dg | doe(e*®/R) exp [in(o — t/a) | <3/|n| (2.12) 


for sufficiently large | n |. It follows that lim,... I(L, r) exists for r < 1 and p > a; in 
fact 


o 


lim I(L,r) = >> r'"'4xI,(T,a)K,(T pe’. (2.13) 
Lo n=—o@ 
If for a particular n, I, is imaginary, then the corresponding term in (2.13) is to be re- 
placed by 22°ir'"' J,,(yn-a)H,” (y,p)e’™ where y, = [k? — (8 — n/a)”]'”. We shall in the 
future write the series as above and leave to the reader the interpretation for imaginary 
ri 
In order to evaluate lim,_,_ lim,... [(L, r), we refer to an earlier paper [4, corollary 
2.1]. It was there shown that for n > 0 





n vl?  —Kiv) on | v .° 
vy E + 4] - oK,,(v) <3 ‘+ n(n—1)] ° (2.14) 
Hence 
Ki(v) n a as n 
K,(v) < v E + "] < v 
and 


K,,(v) < K,,(v9)(v9/v)" for (v > Uo). 
This coupled with the above mentioned fact that 
I,(v)K,(v) < 1/2n 
gives 
I,(T,a)K,(T,p) < (a/p)"/2n (2.15) 
for p > a. Since /_,(T,a)K_,(T,9) = J,(T,@)K,(T,p) the series (2.13) is majorized 


uniformly in 0 < r < 1 by 


Dd’ (a/p)'"'/| 2n |. 
Thus for p > a 
lim lim [(L,r) = >> 4ne*/,(T,a)K,(T,p)- (2.16) 


r-1— Lo n=—@ 
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Finally we shall establish (b), namely, that the lim,.,- Z(L, r) exists uniformly in 
L. To this end, we define 


hr, 0) = P(r, o) for ru, @) Ss ¢, 
=_ for Pe) > ¢, 


ajr,e) = 1 — | do h.(r, c) 


fdr, oc) = h.(r, )/a(r, © 


g(r, o) = [P(r, «) — hr, «))/a(r, ©). 


€ 


It is readily seen that 


| g(r, ¢) do = 1 


0 < f.(r, «) < €/a(r, &) < €/(1 — 2ze) 
and 
a(r, €-) 1 as r—1- 


Now f,(7, «) and g,(r, «) are continuous functions of period 27. Further f,(r, ¢) is 
of total variation 2e/a(r, €) over any interval of length 27. On the other hand g,(r, a) 
is a needle-like function, vanishing outside of an interval of length 6(e, 7) about the origin, 
being otherwise positive valued and enclosing a unit area. It follows that any integral 
of the form J," g.(r, «)w(c) deo differs from w(0) by less than the oscillation of w(c) in 
the non-zero interval of g.(r, «). It follows from known properties of P(r, c) that 6(e, r) 
approaches zero for each easr— 1. 

We write J(L, r) as the sum of two integrals 


oL ‘25 
I(L, r) = a(r, € | dt | da ee" /M) f(r, ¢ — £/a) 
Veak 0 
(2.17) 
oL "25 
+ a(r, ©) | dé | da ee" /Mg lr, ¢ — £/a) 
. L “0 
Expanding f,(r, ¢) in a Fourier series, we obtain 
fr, «) = >. b,(r, ee”, (2.18) 
aaoe 


where 
| b,(r, €-)| < Var (f.)/| m| S 2e/[| n| a(r, 6). 


Substituting this into the first of the integrals of (2.17) and interchanging the order of 
summation and integration (finite limits) gives 


a(r, €) | d¢ | da ee" /R) fr, o — £/a) 
(2.19) 
oL “25 


= D alr, dbatr,) | dg | da e'*(e™*/M) exp fin(o — ¢/a)]. 


n= « Zz /0 











‘0. 3 


r in 
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Finally if we combine this with the inequalities (2.12) and (2.18) we see that the series 
(2.19) is majorized by 


<! 2 
e >’ c/n’ < CL; 


which converges to zero with ¢ uniformly in L for fixed 8 and fixed p > a. 
The second integral in (2.17) approaches [“, dte'*(e"**/R) as r > 1-. Further, we 
note that the oscillation of e“(e"**/§) in the 6(e, r) interval about ¢ = ¢/a is bounded by 


a 


| 
a (e**/R) | Sle, 7) < Cre, r)/R?. 


Hence 


oL Qn L | 
| dt [ da e*(e"* /Rg Ar, ¢ — ¢/a) — [ dg e'**(e'**/R) | 
J-L “0 L ! 


< C,a(,7) | dt = Cr86,7) 


uniformly in L. We have finally that 


| I(L, r)/a(r, &) — lim I(L, r) | < Cye/(1 — 2re) + C3 8(e, r). 
pu] 
Since ¢ is arbitrary and 6(e, r) approaches zero as r — 1 , it follows that lim,_,_ J(Z, r) 
exists uniformly in L. This concludes the proof of (b). Combining this with (a) it follows 
from a well known theorem on interchange of limits that 
IT = lim lim J(L, r) = lim lim I(L, r) 
L 


72 r—1— r>l— Loo 


(2.20) 
= pe A4ne’"* I(T, a)K,(T,,p) 


for p > a. 
It should be remarked that for a point inside of the helix, i.e., p < a, the vector 


potential components can be obtained by simply interchanging a and p since these two 
variables enter symmetrically in the integral J. 

3. Electromagnetic field expressions. In the previous section it was shown that the 
integral J converges for p ~ a. We now define 


Ki, ¢) = [az e*(e**/R), (3.1) 


where R = [¢? + a’ + p” — 2ap cos (¢ — ¢£/a)]'” and 8 is real. It remains to show that 
this integral converges for p = a provided that ¢ ¥ 0 (that is for a point not on the helix 
itself). This will be a by-product of our investigation into the derivatives of [(p, ¢) 
with respect to p and ¢. 

For any pair of values (p, ¢) different from (a, 0) there is a neighborhood of this 
point in which R is bounded away from zero. Further RP is essentially equal to | ¢ | for 


large | ¢|. Hence for the points (p, ¢) of such neighborhood any integral majorized by 
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M/R*, (u > 1), will converge uniformly. This is clearly the case for all derivatives of 
the integrand of I(p, g) with respect to p and ¢. It follows that 


R 


| dp | dé d,|e mi R)] = | dt fe" [(e"** R) — (e" Ro) ]}, (3.2) 


where R, is obtained from FR by replacing p by po . 

Suppose now that pp = a whereas p ¥ a. Then since I[(p, ¢) converges, the right- 
hand side of (3.2) can be written as the sum of two convergent integrals; hence /(a, ¢) 
converges for ¢ ~ 0. Whether pp = a or not, the integrand of the outer integral on the 
left-hand side of (3.2) is continuous. Hence 


\ 


0,1(p, ¢) = 9,[I(p, ¢) — I(po , ¢)] = 9, | dp | dé dle aed R)Ip 


=| dg a,le*(e**/R)). 

It is clear that a similar argument applies to all derivatives of [(p, ¢) with respect 
to p and ¢. Further we now know that J converges for all points off of the wire. This is 
clearly independent of the coordinate system used to express the point. It follows there- 
fore by precisely the same argument as above that the derivatives for the components 
of II with respect to z, y, and z can be obtained by differentiating under the integral 
sign. 

In order to evaluate the derivatives of J(p, ¢) for p ¥ a, we recall the inequality 
(2.15) which shows that the series expansion for J(p, ¢) [see (2.20)] differentiated term- 
wise with respect to yg, is majorized by . (p/a)’” for p < aand by a (a/p)” forp>a. 
It follows that 0,/(p, ¢) can be obtained by termwise differentiation of the series ex- 
pansion. 

A similar result holds for 0,/(p, ¢). Here we make use of the inequality (2.14) for 


F 2 1/2 
es : n v 
K'(v) | < K,(v) — | 1 + ——— ! 
v n(n — 1) 


For v > I,a, the functions (| n |/v)[1 + v?/n(n — 1)]'” are uniformly bounded by some 


the case p > a. Thus 


positive number M. Hence 


< T,M1,(T,a)K,(T,.p) < M,(a/p)" 


l,J,(T,,a)K/(T,,p) 


so that the series expansion for J(p, ¢) can be differentiated termwise with respect to p 


in case p > a. 
Finally p < a, we make use of the corresponding inequality for the J, function 


(n > 0) (see [4, corollary 1.1]) 


n v ae Ae ‘ 

3} 1 7a < - = <3] 1 “3 ; (3.3 

v | + n(n + 5 | —2 Lian | + n ) 
For v < T,a, the functions (n/v)[1 + (v/n)’]'” are uniformly bounded by M’/v for 


some positive number M’. Hence 


l2/(T,.e)K,(T,a) | < (M! ‘p)(p/a)'"', 
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so that again the series expansion can be differentiated termwise with respect to p in 
the case p < a. 

It can now be readily shown that the components of the vector potential (and hence 
the electromagnetic field components) satisfy the wave equation 


A® — €u0,,P = 0 


at all points (z, y, z) not on the helix. In fact, since we can differentiate under the 
integral sign this amounts simply to showing that the quantity e~‘“‘(e“’/r) where 
r= [(x — £)? + (y — n)’ + (2 — §)’]” satisfies the wave equation (Note that k* = 
euw ). 

We shall conclude this section by obtaining explicit expressions for the field com- 
ponents from the formulae (1.2) and (1.6). We first transform the integrals of (1.6) as 
in (1.7). Then we have 


M, = aFe™ | dg e(e™/1) sin ( — z/a — t/a), 


ll, = aFe*” [ dt e'*(e"*"/r) cos (vp — z/a — t/a), (3.4) 


e 


M, =aFe™ | dge*e/n), 


where 
r= [f° + a’ + p” — 2ap cos (y — 2z/a — ¢/a)]'”. 


We shall compute E explicitly. Now 


] 
Vi = 0, + all, + = dally + al, 
8,11, = aFe** | dt e*r"[p — a cos (vp — zya — ¢/a)] sin (y — z/a — £/a)d,(e"7r), 
0.1, = aFe** [ dg e*r-"lap sin (9 — z/a — £/a) cos (y — 2/a — ¢/a)d,(e"7r) 
— sin (y — z/a — f/a)(e"7r)], 
0.01, = isl, + aFe™ [ dt e'*(ar)'[—ap sin (p — z/a — ¢/a)]d,(e"*7r). 
Combining the terms in V -II and making the appropriate cancellations gives 
VII = <6il, (3.5) 


Since E = V(V-Tl) + A’, we have for p > a 
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E, = 180,11, + kl, 


= 4riBaFe'* >, r,/,(T,,a)K,(T,,p) exp [in(g — z/a)] 
(3.6a) 


— i2nk’aFe™ >> [I,-,(T,a)Ky-1(T np) 


n 


— J,.,(T,a)K,.,(T,,p)] exp [in(e — 2/a)], 


= 4riBaF p ‘e'”’ } inl, (T,,a)K,(T,,p) exp [in(g — z/a) 


; (3.6b) 
+ Qnk’aFe'”’ i [7,-:(1,a)K,-.(T,p) 
+ J,.,(T,a)K,.:(1,,p)] exp [in(g — z/a)], 
E, = ipd.Jl, + kl, 
= 4r(k? — B’)aFe'* 4 I,(T,a)K,(T,,) exp [in(g — z/a)] (3.6c) 


. » tBz . " - 1 . 
— 4niBFe ie in I,(T,a)K,(1T,,p) exp [in(g — z/a)]. 

The electric field components for p < a are very similar. One can obtain the magnetic 
field components directly from 


H = —we V X I. 


Since we shall have no occasion to use these in the present paper, we refrain from writing 
them out. 

4. Determination of the propagation constant. In this section we shall determine the 
propagation constant 8 so that the electric field along the helix is ‘“essentially’’ zero. 
Since we assume the helix to be a perfect conductor this is the proper boundary condition 
for our problem. The electric field in the helical direction at any point (p, ¢, 2) is 


E, = (a/[a’ + a’]'”)E, + (a/[a* + a@’)')E, . (4.1) 


In order to approach the helix radially, we set ¢ = z/a in (4.1) and try to determine ~ 
so that lim,.., Z, = 0. We will suppose that p > a; the case p < a is entirely similar. 

Because of the singular nature of the field about an infinitely thin conductor we 
cannot actually choose 6 so that lim,., EH, = 0. Since E becomes infinite as p — a, the 
best we can do is to determine the proper phase velocity, 8) , so that 





lim - — 0. (4.2) 


E,( B | # | Bo ) 


pa 





ic 
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For ¢ = z/a and p > a, EF, can be expressed as 


fa’ +- a |’ “2, = Ama’ (k’? — 8°) Fe” pi I,(T,a)K,(T,.p) 


+ Qra*k?Fe'* > [1,_,(T,a)K,-(0.p) + 1,..(0,a)Ky+.(T,0)] (4.3) 


—4miBa(l — a, p)Fe™ yi in I,(T,a)K, (1, ). 


n=—o 


We shall first prove that the last term in (4.3) approaches zero as p — a. In its 
original integral form this term becomes 


J(p) = (1 — ayp)(ap)~"d,[1(p, ¢)]eas = (1 — a/p) | eR sin ¢ Ag(e"™®/R) dt 


aside from a factor independent of p. Here 


R =([¢ + a + pp? — 2ap cos t/a)” 


wpe 


is an even function of ¢. The e'” can be expanded as 
e'*® = cos Bt + isin fF. 
Since [cos B¢R™ sin ¢ d2(e*”/R)] is odd, its integral vanishes. 
Hence 
J(p) = (1 — a/p) [ R~' sin Bf sin ¢ A,(e""",R) dt. 
We note that 
Is) <i + @— op)? SR. 
It follows that there exist constants c, and c, such that 


d,(e"*/R) | < c,/| ¢ | for Tae 2 


< cn [¢? + (a — p)’ for l¢|<1 and |p-—a|<l. 


Therefore we have 
al 


J(p) < (1 -— ¢ a2 | dé C,/ - + 2 | dt Boot? [¢? .. (a — ar} 


= (1 — a/p){2c, + 2Bc,[log {(1 + [1 + (a — p)’]'”)/(p — a@)} 


— {1+ (@— p)*}~"}}, 
which approaches zero as p — a’. 
In determining the limit of EZ, as p — a*, we can therefore omit the last term in 
(4.3). We denote by E{ the thus abbreviated 2, . We now define the function 


2(8, :) = ¥ 1AP.a)K,(P 0). (4.4) 


n=—@ 
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Then for g = z/a and p > a, we have 


" ie 
Ki = si7a e'*[a*(k? — 8) Q(B, p) 


; la” +a|” (4.5) 
+ (a°k*/2){Q(B — 1/a, p) + 2(8 + 1/a, p)}]. 


It turns out that 
lim Q(B, p) =— Oo, (4.6) 


so that we can hardly expect to obtain a 8 for which lim,.,+ £, = 0. However we shall 


also prove that 


lim Q(B, , p)/Q(B., p) = 1. (4.7) 
Hence 
lim E(B) a’ (k* — B°) +- a’k? 
0 Se ws pe Ea; 
pat L(B,) a(k” — Bi) + ak 
and for 8; ~ 8°. This vanishes if and only if 
a’(k? — B’) + a°k’ = 0. 
In other words the correct propagation constant 8» is simply 
re (4.8) 


Bo = +k{a* + a’) 
We now prove (4.6). In order to do this we make use of the inequalities (2.14) and 
(3.3), from which it follows that 
K/(v) I’(v) n ‘| v | | v | } 
x — ee ee me 8d | 7 . —|1 me : 
vK,,(v) vI,,(v) v \ + n(n — 1) + n(n + 1) 
The expression on the right is monotonic decreasing achieving its maximum of 1/(n" — 1) 
at v = 0. Hence 
(1,(v)K,(v))’ 


x —v/(n? — 1 
I,(v)K,,(v) oa 


and since /,,(0)K,(0) = 1/(2n), we obtain 
(2n)~' > I,(v)K,(v) > (2n) “exp {—v*/[2(n? — 1)]}. 


Finally for sufficiently large | n |, 


7 \2 2)1/2 ! 
r, = [(8 — n/a)” —k])”" >cln| | 

and hence 
I 


I,(1,a)K,(T',a) > | 2n|~" exp (—c’). 


Since 2(8, p) consists of positive continuous terms, monotonically increasing as p — a’, 





we have 


lim (8, p) = 2(B, a) > (1/2) exp (—c’) >» ‘iat, 


n=—2 


p 


which diverges. 
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We shall next establish (4.7) also by means of the inequalities (2.14) and (3.3). 
Now clearly 


O,[1,(va)K,(vp)] _ Pe) o is Ki(vp) 
I,(va)K,(vp) ~ “ I,(va) ' ° K,(vp) 


7 { | J I/(va) | - FE Ki |\ 
= ?\4 | val,(va) oK,,(v) 


It follows from (2.14) and (3.3) that for n > 1 and v > 0 


d,[I,(va)K, vp)] — f, [i a) |” — :| wo" "| 
I,(va)K,(vp) Se oa)? _? Po ry + 


yy\2 | 1/2 yn\2 ]1/2 
enffi+9T'-[+ 9} <0 
v n n 


for p > a. Likewise (n > 1 and v > 0) 


0,{I,,(va)K,,(vp) ‘3 | (va)* - | (vp)? oo |" 
ee > hs - - 
I,(va)K,(vp) — "ly i+ n(n my 1) v ‘+ n(n — 1) 


It is easily shown that the function within the brackets is negative monotonic increasing 











for p > a, achieving a minimum of 


(. n(p- yt +a) ” joe 


—] 
Hence 
a,[I.(0a)K,(v9)] . _ 
T,(va)K,(vp)  — 
where 





Wy _ 


nip —a)+(p +2) 
n — 1 
If we now integrate the above logarithmic derivative from.v > 0 tov + A > 0 we obtain 


forn > 1 








: Wn A . I,[(v + A)a]K,[(v + A)p] 
exp { 9 [v + A) —v I) < <2 Pl<1 for A>0O (4.9a) 
and 
I,[@ + A)a]K,[v + Ae] - {.. in a ")) 
1< T.(va)K.(vp) < exp 9 [wv + A) v'| for A <0. (4.9b) 


Finally because J_,(va)K_,(vp) = I,(va)K,(vp), the above inequalities hold forn < —1 
if we replace n in w, by | n |. We can combine (4.9a) and (4.9b) into 


(pan) mae [(v + A)a]K,[(v + A)p] ,€(p.m) 
é < T.(va)K (op) <e (4.10) 
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for all|n| > 1,v > 0, (v + A) > O, and p > a, where 


in | (p- — @) +p + a’ | (y + A)? - y* | 


(p, = 
e(p, n) = oe | 


It will be seen from the expression for (8, p) that the appropriate values for v and 


(v + A) in I,(va)K,(vp) are 


> | (2. = ny ~ «| = cag 
a 


and 


Hence 
(v + A)? — v2 = Br — Bi — 2n(B. — B,)/a. 


It is clear that ¢(p, n) is the sum of two parts one of which is of the order of (p — a) 
and the other of the order of 1/(| n |). Hence given 8, , 6, , and » > O, there exists an 
N and a6 > O such that for|n| > NandO <p-—a<6 


e(p, n) < 7». 


The finite sums be I,(1,’a)K,(T\"’p), remain uniformly bounded as p — a’. Hence 


coe Q(8; p) = N * . a 
Q(B. , p) Q(B. , p) 
N f N f | 2 
| 14 sf 
SU Ett es 
Q(B. , p) Q2(B2 , p) / 
where { } stands for the running term 
(1.(r aK, (Te p) — 1,(P."a)K, (0, p)}3 


thus 


QB; , p) 


N ae a! —— — a. ee 
Q(B. ,p)| ~ (8, p) +t Q(B. , p) 
It follows that 
— Q(B, , p) 
lim j|}1—-— <O0+e’-—e"’ 
p-a* Q(B. ’ p) ; 


and since 7 is arbitrarily small, we have 
lim Q(8, , 9) a] 
pat (8, ’ p) 


5. Conclusion. The net result of our investigation has been to obtain expressions for 
the electromagnetic field which satisfy the Maxwell equations away from the helix and 
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which satisfy the appropriate boundary conditions along the helix. The propagation 
constant in the axial direction was found (4.8) to be 


By = +kfa? + a’]'7a7 (5.1) 


which corresponds to precisely the phase velocity one would obtain if the wave traveled 
along the helix itself with the free space velocity of propagation. 
In case 


(8, — n/a)’ > k’ (5.2) 


for all n (n = 0, +1, +2, ---) none of the terms in the expansion of the field equations 
are unmodified Hankel functions. As a consequence there is no net radial flow of energy. 
On the other hand for a given helix there is a series of frequency bands for which the 
inequality (5.2) is not satisfied for all n; in this case a finite number of Hankel function 
terms do appear in the field expressions. It can be shown that these terms do cause a 
net radial flow of energy. One might expect this phenomena to be reflected in the propa- 
gation constant 8) ; that is, one would expect 8, to be complex. More precisely, unless 
there is a source of energy within the wire itself, one would expect 


| Energy radiated per unit distance in the z-direction 


UBo) = aie — _ : : : 
Im(Bo) 2 Energy flow in the z-direction 


(5.3) 





However the computed value of 6, (5.1) is actually consistent with (5.3) for it can be 
shown that whereas the energy radiated per unit distance in the z-direction is finite, 
the energy flow in the z-direction is infinite. As a matter of fact not only is the total 
contribution of the set of terms satisfying (5.2) infinite, but the contribution of each of 
the radiating terms to the energy flow in the z-direction is itself infinite. 

In spite of the physical reasonableness of our result, there is, in the case of the radi- 
ating helix, reason to be dissatisfied with it. If (5.2) is satisfied for all n, the expansion 
which we have obtained can be used for helices with conductors of finite thickness. One 
could, for instance, replace the delta function of section 2 by a finite needle function 
(say P(r) , «) with fixed r, < 1). The helix would then be constructed from infinitely 
thin tape. For instance II, would be of the type (p > a) 


w 


lI, = aFe’* > a,I,(T,a)K,(T,) exp [in(g — z/a)] 
where the (a,) are Fourier coefficients of the needle function. Repeating the above de- 
velopment one would obtain a real 8, for the a,’s sufficiently close to 1. However, in 
the radiating case, one would not be able to obtain a solution in this way. In this case, 
no real 8 would allow the boundary conditions to be satisfied whereas the entire de- 
velopment (starting with the vector potential integrals) is meaningless for complex 8. 
This is also reflected by the fact that a complex 8 which satisfies the boundary conditions 
would make the arguments of all of the Hankel functions complex so that the field 
would decay exponentially. As a consequence we would have no radiating terms which 
would be contrary to (5.3). Thus in case (5.2) is not satisfied by all n, one would expect 
a real helical wire to radiate. However one would not expect the field expressions for 
finite conductors to have expansions of the type which we have obtained for this case. 








246 R. S. PHILLIPS [Vol. VIII, No. 3 


REFERENCES 


1. Franz Ollendorff, Die Grundlagen der Hochfrequenztechnik, Springer, Berlin, 1926, pp. 79-87. 

2. J. R. Pierce, Theory of the beam-type traveling-wave tube, Proceedings of the Institute of Radio 
Engineers 35, 111-123 (1947). 

3. R. 8S. Phillips and Henry Malin, A helical wave guide, New York University Mathematical Re- 


search Group, Report No. 170-3, Army Air Forces, Watson Laboratories, 1947. 
4, R. 8S. Phillips and Henry Malin, Bessel function approximations, to be published in the Amer. J. 


Math. 
5. S. Ramo and J. R. Whinnery, Fields and waves in modern radio, Wiley and Sons, New York, 1944. 


6. G. N. Watson, Theory of Bessel functions, MacMillan, New York, 1944. 
7. A. Zygmund, Trigonometrical series, Monografje Matematyczne, Warsaw, 1935. 











247 


TRANSIENTS IN MULTIPLY PERIODIC NON-LINEAR SYSTEMS* 


BY 
FRANK E. BOTHWELL 
The University of Chicago 


1. Introduction. The transient behavior of the general dynamical system of n degrees 
of freedom, subject to certain restrictions, will be treated. The system must be (a) nearly 
linear, (b) nearly conservative, (c) nearly autonomous; the deviation from the case of 
complete linearity, conservativeness, and autonomy being indicated by the magnitude 
of a parameter yz. It is further assumed that the explicit dependence of the non-auton- 
omous terms on time is in the form of sinusoids of frequencies Q, , 2, , «++ , Qn . 

The general system is transformed to normal coordinates and thence to elliptical co- 
ordinates, the resulting form being convenient for the development of the solution for 
small values of uw. An auxiliary set of differential equations is so constructed that its 
singular points determine the limit cycles of the original system for the limiting case, 
u — 0. The stability of each limit cycle in the original system is the same as the stability 
of the corresponding singular point in the auxiliary system. Furthermore, the coordinates 
of the auxiliary set are the amplitudes of the n fundamental frequency components of 
the original system, so that the motion of the representative point in the phase space 
of the auxiliary system describes approximately the transient behavior of the amplitudes 
of the fundamental frequency components of the original system for small values of the 
parameter yz. 

The general cases of internal and external resonances, in which two or more funda- 
mental frequency components are commensurate, are treated. Finally, the method is 
applied to a problem of practical interest. 

2. The method. Under restrictions (a), (b), and (c), the equations which describe the 


general system are 


pm (a, , — bus) + ug; = 0, ¢=1,2,--- ,a. (1) 
i=1 
The quantities g; are functions of the parameter u, the independent variable t, and the 
n dependent variables y; and their various derivatives. The final restriction on the 
system (1) is that there exist a closed region D, including the origin, in which each g; is 
bounded for all values of ¢ and for sufficiently small values of x. 

By means of a suitable linear transformation [1] 


es 2, CiBe 5 7 = l, 2, ee aa 


k=1 


Eqs. (1) may be transformed to normal coordinates. Thus 


d’x, 2 ‘ 
“ae + et + uf; = 0, ; = 1.2, ** 5 Me (2) 
The constants w, are the roots of the characteristic equation of the linear system which 


*Received July 22, 1949. 
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occurs when uy is set equal to zero in (1). Since the kinetic energy of the linear system is 
a positive definite form, the constants w; are positive [1]. 

The classical approach to equations (2) involves the expansion of the normal modes 
x; and their fundamental frequencies in power series in u. This procedure is an extension 
of the method presented by Kryloff and Bogoliuboff for a system of a single degree of 
freedom [2]. Concerning this extension, only one remark will be made. In the system of 
a single degree of freedom, secular terms in the series solution may be avoided by satisfy- 
ing two equations at each stage of the process. One of the two equations determines the 
amplitude of the corresponding term in the series, and the other equation determines 
the contribution of the corresponding term to the fundamental frequency of oscillation. 
In the case of a multiply periodic system, each single equation is replaced by a set of n 
simultaneous non-linear equations. 

For the purpose of this note, it is convenient to use elliptical coordinates defined by 


the equations 


ZX; = 7; cos 0; 
7=1,2,--- ,n. (3) 
dx; —* 
- = q,7r,; sin 6; 
dt 
In these coordinates, Eqs. (2) become 
dr mm ; 
—+ = —— f; sin 0, 
dt w; * 
2 = I, 2, - ,n (4) 
dé m 
— = —w; — f; cos 8; 
dt m w.7; * ' 
By means of the change of variable, 
T = ut, 
Eqs. (4) may be written 
dr ar 
— = —— f,sin 6; 
d7 a“ 
= 1,2, +o 2”. (5) 
dé, w; = 
—_ = —— — — f; cos 6; 
dl m w,7; 


Inasmuch as each g, is bounded in D, each f; must be bounded in D. Since the quanti- 
ties f; are functions of the variables x; and their derivatives, and since the latter are 
periodic functions of the variables 6; , the right sides of (5) are periodic functions of 
each variable 6; . Furthermore, by hypothesis, each f; is a periodic function of each 
variable @; where 


° 
d= -—T+ 0, t= 1,2, +++ ,m. 
bb 


As the parameter yu is decreased, the right sides of (5) take the form of certain slowly 
varying functions of 7 (due to the variation of r; with time) modulated by rapidly 
varying quasi-periodic functions of 7'(due to the variations of 0; and ¢,; with time). 
Furthermore, the quasi-periods of the modulation approach zero with u. 





ly 
ly 
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Let K;(r, , 2, +++ T, , #) represent the unmodulated portion of f; sin 6; , defined by 
» 
K fe) = tie. : itil C0t oe (6) 
Pox “0 


in which the variables r; are considered constant in the integration. When evaluating 
the integral (6) for the limiting case u — 0, the integrand is a function of each r,; (con- 
sidered constant in the integration), a sinusoidal function of each 6; (equal to —w,t + Oi 
in the limiting case), and a sinusoidal function of each ¢; (equal to —Q,t + 4). 

Now the quantities 


G, = f, sin 6; — K,(u), += 1,2, ---,n, 


are quasi-periodic functions of 7’, the quasi-periods of which approach zero with uz. 
Furthermore, each G; is bounded in D, so that the limit as u — 0 of the indefinite integral 
of G; is identically zero. Therefore, for the limiting case, nu — 0, Eqs. (5) may be written 


dr, - 
= —-—— += 1,2,---,n 
d7 ox” Dae is 


and for small values of yu, (4) may be approximated by 


dr; ae . on 
> las Kw > ¢= 1,2, --- ,n. (7) 
The system of n simultaneous non-linear differential equations (7) comprises the set 
which is auxiliary to the system (4). In fact, the variables of (7) are the amplitudes of 
the set (4). Therefore, the singular points of the former, defined by the equations 
K;(rio »T205 *** » Tro ’ 0) = 0, t = 1, 2, iia (8) 
fix the amplitudes of oscillation on the limit cycles of (4). Furthermore, the stability of 
each singular point of (7) is the same as the stability of the corresponding limit cycle 
of (4). 
Let us consider first the non-resonant case, that is, the case in which no two of the 
set of frequencies w; and 2; are commensurate. In this case, (6) may be replaced by the 


expressions 
; | a2 a2 : 
KO =——< | | - ] fisin 6d, d0, --- d0, dg, dbz -+- dbn, — 9) 
(2r)""" 40 0 Jo 


¢=1,2,--- ,n. 
Thus the values of AK; , the singular points of the auxiliary system (7), and hence the 
limit cycles of (4) for small u are independent of the phases @;) and ¢,o of the autoperiodic 
and heteroperiodic oscillations. 
3. Internal and external resonances. If two frequences, say w, and w, are commen- 
surate, that is, if there exist two positive integers n, and n, such that 


NW, = NoW. , 


the system will be said to have a simple resonance of order n,/m2 or n2/n, , Whichever 
is smaller. On the other hand, if there exist k positive integers, n, , m2, +++ , m , Such that 


NW, = NeW. = °° = UA, (10) 
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the system will be said to have a complex resonance of degree k, there being associated 
with the resonance k — 1 orders, n2/n, , %3/n, , -** , M%/n, , Where n, is the largest of 
the integers. 

In this case, the integral (6) may be replaced by 


1 a2x a2r a2r 


K,(0) = ari | | | f, sin 0; dx d0,., +++ d0, db, +++ dbn , (11) 


Zn “0 “0 “0 


a= 1, 2, . i 

in which 

N 

6,:=>—x2+ 600, t=1,2,°-:,k, 

n, 

where 
N = I] N; 
1=1 

Thus the integration on x from zero to 27 insures that 6, , 0. , --- 6, range over integral 


numbers of their respective periods. 

If only frequencies w; and not Q; are involved in (10), the resonance is said to be 
internal. If only frequencies 2; and not w,; are involved, the resonance will be called 
external. If both w; and Q; are involved, the resonance will be called compound. 

In the general case, one may have a system with several complex compound reso- 
nances. In any case, it is always possible to find a set of integers (sufficiently large) 
such that (10) is approximately true. Therefore, the resonance case makes sense only 
provided the integers are all small. The values of the integrals (11) in many cases depend 
upon the phases @;5 and ¢io . 

4, Example. Consider the circuit of Fig. 1, containing the non-linear resistor 


R= —p(l —yT). (12) 


The resistance R is composed of a passive resistor in series with the negative resistance 
of a vacuum tube. The passive resistor is so chosen that the total resistance of the circuit 
at zero current (that is, —) is small and negative. The term y’J’ is due to saturation 
of the tube. 


R L, L, 


Fs i) - (1 a 


Fic. 1. 
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The two degrees of freedom of the circuit of Fig. 1 may be uncoupled by the trans- 


formation to normal coordinates 


Rh=K+1; 


I, = al; + B12, 
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in which 


1 
a fiieaaed Qu [w32 = wn + { (wie ro wi)” + 42, }'/*], 
12 
1 
B= Qs, we. apd wr a { (wee a wi)” + 4erioe9,}'”*], 
12 


where 


| i 


mG | 1 
of = P(E +4), wis = ? 


SS 
° 


ae 
%! 


n= 2(2 1) _— 
®o22 = : a. a Oa = 


The resulting circuit is shown in Fig. 2. 
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The equivalent circuit elements are given by the expressions 





tats ———— = 

L="S4h,, C= 77s 
(13) 

eet. ee B  — - @ 

Li = e L, ; C3 — Le 


where 


2 1 9 9 9 9 9 9 9 D 
OF = 5 let + wr + {Gi — we)’ + deine }'”*], 


9 1 9 9 9 ” 2 2 2 6 
03 = 5 [wir + wi» — {( — wx)” + Aes oo1}'/? J. 


In the circuit of Fig. 2, the inductances and capacitances are positive according to 
Eqs. (13). 
It is convenient to make the substitution 


z,=7;, += 1, 2. 
The Kirchhoff equations for Fig. 2 are 


d "i 


dx, ds) 
dt’ 


+ Ox; + uf (x, ’ Xe Eee 


may t= (14) 
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in which 


Lif. = (er + ay? — 1( + aa) = Lif. . (15) 
: dt dt 
Equations (14) together with (15) are an extension of Van der Pol’s equation. 
For the case of u very small, one may proceed directly to equations (7) after evaluating 
the integrals (6). If Q, and Q, are not commensurate (non-resonant case), the corre- 
sponding integrals (9) are 


Ao 1 a2r a2e ; : 
K {G@) = er | | [(r; cos 0, + r. cos 62)) — 1](w,7, sin 8, 
®L:Jo Jo (16) 
+ wr. sin 6.) sin 6; dé, dé, , i=1,2 
so that equations (7) take the form 

dr M 2 2 

rs = 377 r(4 — rj — 2r2) 
(17) 

dr he 


dt ~ 8L r(4 — rz; — 2ri). 

The singular points of (17) occur at (0,0), (2,0), (0,2), and (2/3'”, 2/3’). They 
are respectively an unstable star point, a stable nodal point, a stable nodal point, and 
a saddle point. The corresponding limit cycles of (14) have the sets of amplitudes (0,0), 
2,0), (0,2), and (2/3'””, 2/3'*) and are respectively unstable, stable, stable, and un- 
stable. 

The transient behavior of the amplitudes of the normal oscillations of (14) can be 
shown in the 7,7, phase plane of the system (17). The differential equation of the phase 
plane trajectories is obtained by dividing the second Eq. (17) by the first. Thus 


= we 9 Tee : 
dr, [7,4 —-—17r; —2 


t 


dr. Lir, 4 — r> — 2r 
>, 


to 


For the case Li = Li , the solution is 


(rr; — ri)’ = iE _ rt +- riri(c — 2lIn ") | (18) 
t 


c being the constant of integration. The trajectories are shown in Fig. 3. 
If 2, and Q, are commensurate, that is, if there exist two integers, n, and n, , such 
that 


n,Q, = nL, , 


> 


equation (11) must be used to evaluate K, and K, . The integrals are easily evaluated 
and yield four cases depending upon the ratio n/n, . If the ratio is unity, the auxiliary 
set of differential Eqs. (7) contains the synchronous phase 6,5 — 620 . If the ratio is 1/2 
or 3, the auxiliary set contains the synchronous phase 36,5 — 029 or 302. — 0. . Any 
other ratio yields equations identical with Eqs. (17) describing the non-resonant case. 
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No. 3 
The location and stability of the limit cycles of the original system, as well as the 
transient behavior of the amplitudes of the normal modes, are strongly affected by the 
(15) order of resonance. In the above example for the case of internal resonance of order 
5 
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ae ars ry, a 
yrre- ae. 
[p 


(16) LS \ Ps 


} PA 
| WZ RN 
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1 
Fig. 3. Transient variation of the amplitudes of the normal modes of a doubly periodic system. 
ich 
one-third (n, = 3n,), there are four limit cycles. Three of the limit cycles have the re- 
spective amplitudes (0,0), (2,0), and (0,2), and correspond respectively to an unstable 
star point, a stable nodal point, and a stable nodal point in the auxiliary system. The 
ied remaining singular point is a saddle point; its location is a function of the synchronous 
ry phase 
9 
i 9 = 3610 — 420, 
ny 


and is shown in Figure 4. 
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Except for internal resonance of order one or one-third, the auxiliary equations are 
the same as for the non-resonant case, and hence the transient motion is not affected. 
The effect of resonance on the location and stability of limit cycles is complicated even 
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Fic. 4. Effect of synchronous phase, % , on the amplitudes of oscillation of the normal modes of a 
doubly periodic system. 

in this simple case of a doubly periodic system, and will not be discussed in more detail 

here. 
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PSEUDO CLOSED TRAJECTORIES IN THE FAMILY OF TRAJECTORIES 
DEFINED BY A SYSTEM OF DIFFERENTIAL EQUATIONS* 


BY 
L. A. MacCOLL 
Bell Telephone Laboratories 


1. Introduction. This paper is concerned with certain simple closed curves, here called 
pseudo closed trajectories, which play an important part in determining the topological 
properties of the family of trajectories (or characteristics) defined by a system of differ- 
ential equations of the form 


f= Xz), = Ve,y. (1) 
Some of these curves are considered in a rather incidental way in the writings of Poin- 
caré.** However, the full concept of pseudo closed trajectories does not seem to have 
been discussed explicitly heretofore. 

We assume that all of the variables in the equations (1) are real, that the functions 
X and Y are continuous in an open connected region F in the zy-plane, that these func- 
tions satisfy Lipschitz conditions locally throughout R, and that the real curves X = 0 
and Y = 0, if they exist, have only simple intersections. We assume that these points 
of intersection have no point of condensation in the region R. We also assume that if 
the functions X and Y are not linear, they satisfy conditions which are sufficient to 
insure that the topological properties of the family of trajectories in a small neighborhood 
of any singular point (point of intersection of the curves X = 0, Y = 0) other than a 
center are the same as in the case obtained by replacing X and Y by their linear ap- 
proximations. The last assumption is satisfied if X and Y are, for instance, of class C° 
in R. 

In brief, our assumptions are simply the ones that are usually employed in discussions 
of the family of trajectories. The necessary information about the implications of the 
assumptions is readily available in the literature. 

Since the notion of a trajectory is a familiar one, we have used the term so far without 
any explanations. However, in order to avoid the danger of future ambiguities, it will 
be well now to define the term explicitly. 

Let x = o(t), y = V(t) be a solution of the system of equations (1). It is understood 
that the solution is defined over a certain maximum (finite or infinite) open interval 
I‘® on the t-axis. 

As t describes J‘” once in the sense of increasing t, the point (x, y) = (¢(é), ¥() 
describes a certain curve 7’ in R. We call T a trajectory. 


*Received Aug. 3, 1949. 

**See, for instance, Oeuvres de Henri Poincaré, t. I, p. 56 et seq. Poincaré’s considerations concerning 
these curves are limited to pseudo closed trajectories which have only exterior vertices or only interior 
vertices, and to the possibility of trajectories being asymptotic to such pseudo closed trajectories. 
Exterior and interior vertices are defined in the next section. 

+A good discussion, covering most of the points of interest, is given in the recent book by Tricomi, 
Equazioni differenziali, Einaudi, Torino, 1948. 
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Any singular point is itself a trajectory, a point trajectory; and these are the only 
point trajectories. A trajectory 7’ which is not a point trajectory is oriented, the positive 
sense on T' being that in which the curve is described as ¢ increases. 

2. Pseudo closed trajectories. Under our assumptions we may have just four kinds of 
singular points: nodes, foci, centers and saddle points. Since the topological properties 
of the family of trajectories in the neighborhood of any one of these singular points are 
well known, it is unnecessary to review them extensively here. 

We recall that a node or focus is said to be stable or unstable according as the oriented 
ares of trajectories in the neighborhood proceed toward or away from the singular point. 

In the neighborhood of a saddle point there are just two ares of trajectories which 
proceed toward the saddle point. These arcs, considered for the moment as being un- 
oriented, together with the saddle point itself, form a smooth are. Also, there are just 
two ares of trajectories which proceed away from the saddle point; and these, together 
with the saddle point, form a smooth are. The two smooth ares described have distinct 
tangents at the saddle point. These phenomena are illustrated in the typical Fig. 1. 


D Cc 
p 
B 
A 
Fia. 1. 


Elementary examples show that a trajectory may proceed from a saddle point to 
the same or a different saddle point.* We now assume that we have a case in which 
certain trajectories of this kind exist. Let the symbols 7, , T, , --- denote these tra- 
jectories. 

Using a finite set of the trajectories 7; , and also using saddle points as necessary to 
connect these trajectories, we try to construct a closed curve I satisfying the following 
conditions: 


C8 I is oriented, the orientation being that given by the orientations of the 
trajectories 7’; contained in I. 
Cy: I is a simple closed curve. 


Thus, referring to Fig. 1, C, implies that if AP is an arc of T, the succeeding arc is PB 
or PD. Also, if AP is an are of I’, C, implies that the are CP does not belong to I. 

We shall call a curve I satisfying the above conditions a pseudo closed trajectory, 
or briefly a p.c.t. There are families of trajectories which contain no p.c.t. On the other 
hand, familiar or easily constructed examples show that there are cases in which p.c.t. 
do exist. In all of the following we assume that we have a case of the latter kind. 


*See, for instance, Andronow and Chaikin, Theory of oscillations, Princeton, 1949, p. 71. 


iC 
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We observe that different p.c.t. may have isolated points, or entire arcs, in common. 

We shall call the saddle points which lie on a p.c.t. the vertices of the p.c.t. It is 
necessary to draw a distinction between two kind of vertices, which we call exterior 
vertices and interior vertices, respectively. 

Let I be a p.c.t., and P one of its vertices. As a Jordan curve, [ divides its point 
set complement in the plane into a exterior region and an interior region. In the neigh- 
borhood of P we have an are of trajectory proceeding toward P, and an arc of trajectory 
proceeding from P, neither of these arcs belonging to I’. If both of these ares of trajectories 
are in the exterior region, we say that P is an exterior vertex; and if both of these arcs 
are in the interior region, we say that P is an interior vertex.* In the case of the p.c.t. 
lr shown in Fig. 2, P, and P, are interior vertices, and P; is an exterior vertex. 

3. The fundamental theorem concerning pseudo closed trajectories. 

THEOREM. Let T be a p.c.t. such that all points of the interior region determined by T 
belong to R. Let N, , N,, N.,N, , respectively, denote the numbers of nodes, foci, centers, 
and saddle points in the interior region, and let N;, denote the number of interior vertices 
of T. Then we have the relation 


N.+Ne, +N. —-N,=14+N,,. 


The proof of this theorem depends upon the following two familiar theorems. 

T, : Let T be a simple closed plane curve consisting of a finite set of smooth arcs 
meeting at corners. As the running point describes I once in the counterclockwise sense, 
the tangent rotates through + 27 radians, it being understood that at an exterior 
(interior) corner the tangent rotates counterclockwise (clockwise) through the ap- 
propriate exterior (interior) angle. fT 

T, : The Poincaré index** of a simple closed curve which does not pass through any 
singular point, and which is such that the interior region determined by the curve belongs 
to R, is equal to the sum of the numbers of nodes, foci, and centers in the interior region, 
diminished by the number of saddle points in the interior region. 

Let T be a p.c.t. with interior vertices P, , --- , P, . It is understood that these points 
are ranged in the order in which they are encountered when I is described in the sense 
of its natural orientation. 

We choose points (not vertices) Q{ , Qi’, --- , Qf , Qi’ on T so that the typical arc 
Q'Q’ of T contains P; , and no other vertex. At each of the points Q{ , Q/’ we draw the 
normal to I. 

Now we construct a simple closed curve I* as follows. 

(1) For each 7, we join a point on the exterior normal at Q/ to a point on the exterior 
normal at Q/’ by means of an arc of trajectory. The classical continuity theorems, together 
with the known behavior of the trajectories in the neighborhood of a saddle point, insure 
that this construction is possible, provided that the initial point of the are of trajectory 
is taken sufficiently near to I. 


*It is clear that it is impossible for one of these arcs to be in the exterior region and the other in the 
interior region. 

+A rigorous modern proof of this intuitively evident theorem is given by J. Lifshitz, Boletin de la 
Sociedad Matematica Mexicana, 3, 21-25 (1946). 

**Following the usual modern practice, we reckon the indices as the negatives of the indices actually 


defined by Poincaré. 
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(2) Fori = 1, --- , k — 1, we join a point on the interior normal at Q?’ to a point 
on the interior normal at Q{,, by means of an are of trajectory. If the initial point of 
the are is taken sufficiently near I’, this construction is possible, for the reasons mentioned 
above. Similarly, we join a point on the interior normal at Qj’ to a point on the interior 
normal at Q{ by an arc of trajectory. 

(3) We complete I* by joining the end points of the arcs of trajectories that we have 
drawn by short segments of the normals to I at the points Q/ , Q?’. 

For the sake of greater clarity, we emphasize the fact that the curve I'* is supposed 
to be constructed so that each arc of trajectory belonging to I™ lies very close to the 
corresponding are of T (i.e. an are Q{Q}’, an are Q?’Q/,, , or the are Q;’Q{) throughout 
its length. This is possible, because we have constructed the curve so as to avoid the 
ares of trajectories which proceed to and from the vertices of [', and which do not belong 
to I. 

Fig. 2 illustrates the construction of I'* in a simple case. 





Fia. 2. 


The curve I* is closed, by its construction. A little consideration suffices to show 
that if the segments of normals to I’ belonging to I'* are taken sufficiently short, I'* is 
a simple closed curve, and that the interior region determined by I* belongs to R. Also, 
under the same condition, ['* does not pass through any singular point; and the singular 
points contained in the interior region determined by I* are the singular points con- 
tained in the interior region determined by I, plus the interior vertices of T. Therefore, 
by T, above, the index of [* is VN, + N, + VN. — N, — Ni, . 

We now proceed to evaluate the index of I* in a different way. 

The index of I'* is the net number of revolutions performed by the field vector (X, Y) 
in the positive sense when the point (z, y) describes T* once counterclockwise. The 
continuity of the functions X, Y implies that the contribution to the index, conceived 
in this way, coming from the short segments of normals to I is small in absolute value. 
Now consider the contribution coming from the arcs of trajectories belonging to I*. 

We note that at a point of I’, other than a vertex, the field vector does not vanish, 





it 
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and has the direction of the tangent to I. From this, and the continuity of X and Y, 
it follows at once that the contribution to the index of '* coming from the arcs of tra- 
jectories belonging to I'* is nearly equal to the net number of revolutions performed by 
the tangent to T when the point of tangency describes I once counterclockwise. By T’, 
above, this number of revolutions is +1. 

The conclusions of the preceding two paragraphs imply that the index of I'* cannot 
be far different from +1. Hence, since the index is necessarily an integer, the index is 
+1 exactly. 

When we equate the two expressions for the index that we have obtained, we have 
the conclusion of the theorem. 

The preceding reasoning is based upon the hypothesis that T has at least two interior 
vertices. However, the same proof, with some slight and obvious modifications of the 
expository details, applies also when T has only one interior vertex. 

If f has no interior vertices, we draw the normal to IT at a point Q which is not a 
vertex. Then we construct a simple closed curve I'* out of a short segment of the interior 
normal and an are of trajectory which joins the end points of the segment, and which 
lies very close to T throughout its length. Having constructed I'* in this way, we obtain 
the conclusion of the theorem by means of reasoning which is essentially similar to that 
used above. 

4. Study of some simple cases. In the case of a p.c.t. that has no interior vertices, we 
have the relation 


N,+tN,+N. -—N, = 1. 


This is the same as the relation that we have in the case of a closed trajectory (in the 
ordinary strict sense). Consequently, the existence of a p.c.t. without interior vertices 
leads to essentially the same conclusions, in regards to the topological properties of the 
family of trajectories in the interior region, that the existence of a closed trajectory does. 
Since these conclusions have been thoroughly discussed in the literature, we shall not 
consider them here. Instead, we proceed to consider some simple cases in which we have 
p.c.t. with interior vertices. 


CAsE 1. 


Suppose that T is a p.c.t. with just one interior vertex P, and suppose that there 
are no centers, saddle points, or closed trajectories in the interior region determined by 
I’. (It is understood that the interior region belongs to R.) 

The fundamental theorem implies that there are just two singular points in the 
interior region. Either of these points, which we shall denote by the symbols S, , 8; , 
may be a node or a focus. We shall see that when the types of these singular points are 
determined, there remain just two topological possibilities as to the structure of the 
family of trajectories in the interior region. 

In the interior region there is one arc of trajectory which proceeds from P, and one 
are of trajectory which proceeds to P. In the first subcase, illustrated by Fig. 3a, these 
ares belong to the same trajectory 7’. The trajectory 7 and the saddle point P form a 
p.c.t. I’ with just one vertex, which is exterior. The interior region determined by I’ is 
contained in the interior region determined by I’. The fundamental theorem implies 
that one of the points S, , S, , say the point S, , lies in the interior region determined 
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by I’. Then S, must lie in the complement of this region with respect to the interior 
region determined by I. 

Suppose that the ares of trajectories which proceed from and to P, respectively, do 
not belong to the same trajectory. Then the arc of trajectory proceeding from P cannot 
proceed to any singular point on I, for all of the ares of trajectories proceeding to such 
singular points are otherwise accounted for. The are of trajectory cannot approach 
asymptotically to a closed trajectory, for there are no such trajectories in the interior 
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region. Simple considerations, which we omit for the sake of brevity, show that the arc 
of trajectory cannot approach asymptotically to T. It follows, therefore, that the are of 
trajectory must proceed to one of the points S, , S, , say the point S, . Similar reasoning 
shows that the are of trajectory which proceeds to P must proceed from S, . 

We have the situation illustrated in Fig. 3b, which is drawn on the basis of the 
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assumption that S, is a node and S, a focus. The appropriate modifications of the figure 
for the cases in which S, and S, are of other types are obvious. 

We note that in this second subcase S, and S, must be stable and unstable, respec- 
tively. 
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CASE 2. 


Suppose that we have a p.c.t. I! with just two interior vertices, P, and P, . Suppose 
also that the interior region determined by I’, which is understood to belong to R, con- 
tains no centers, saddle points, or closed trajectories. Then the fundamental theorem 
implies that there are just three singular points, S, , S, , and S; , in the interior region. 
Any one of these may be either a node or a focus. For the sake of simplicity we draw 
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the figures for the case in which all three of these singular points are nodes. The modi- 
fications of the figures for the other cases are obvious. We shall see that when the types 
of S, , S, , S; are determined there remain six different possibilities as to the topological 
properties of the family of trajectories in the interior region. 

1. Consider the two ares of trajectories, in the interior region, which proceed from 
P, and P, , respectively. It may be that the are proceeding from P,(P,) and the are of 
trajectory proceeding to P,(P,) belong to the same trajectory 7,(T,). This is the situa- 
tion illustrated in Fig. 4a. Reasoning which is similar to that employed in the preceding 
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case shows that if the subscripts are assigned suitably, the singular points S, , S, , S; 
must be located as shown in the figure. 
2. It may be that the are of trajectory proceeding from P, and the arc of trajectory 


ny 


proceeding to P, belong to the same trajectory 7, , but that the like is not true of the 








262 L. A. MacCOLL [Vol. VIII, No. 3 
ares of trajectories proceeding from and to P, . Then one of the singular points S; , say 
the point S, , must be located in the interior region determined by the p.c.t. I’ formed 
by 7, and P, . Also, we see that the arc of trajectory proceeding from P, must proceed 
to S, or S; say to S;, , and that the arc of trajectory proceeding to P, must proceed 
from S; . This situation is illustrated in Fig. 4b. 

3. It may be that the are of trajectory which proceeds from P, proceeds to P, , 
and that the are of trajectory which proceeds from P, proceeds to P, . This situation is 
illustrated in Fig. 4c. The fundamental theorem implies that the points S, , S, , Sz , 
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with a suitable assignment of the subscripts, must be located with respect to the tra- 
jectories as shown. 

4. It may be that the arc of trajectory which proceeds from P, proceeds to one of 
the singular points S; , say the point S, , and that the arc of trajectory which proceeds 
from P, proceeds to P, . Then it is easily seen that the arc of trajectory which proceeds 
to P, must proceed from another singular point, say S, . The third singular point must 
be located as shown in the illustrative Fig. 4d. 
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5. In each of the four preceding subcases we have at least one are of traectory which 
proceeds from one of the vertices P, , P, to the same or the other vertex. Now suppose 
that this condition is not satisfied. Then the ares of trajectories which proceed from P, 
and P, must proceed to the same or different ones of the points S; . Suppose that the 
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. two ares proceed to the same singular point S, . It is clear that the arcs of trajectories 
proceeding to P, and P, , respectively, must lie on opposite sides of the arc P,S;P, . 
Consequently, these arcs of trajectories must proceed from singular points S, and 8, , 
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= which are located as shown in the illustrative Fig. 4e. 
° ~s 
iS 4 \ 
/ s 
! \ 
: I 
t S, 
3; / 
% f° 
\ 7 
* 2 
~ = 
= 
™ ~ee ae = - 
Fia. 4e. 
6. Suppose, on the contrary, that the ares of trajectories proceeding from P, and 
P, proceed to S, and S, , respectively. Then it is easily seen that the arcs of trajectories 
proceeding to P,; and P, must proceed from S, . This subcase gives us a variant of Fig. 4e 
which we do not need to exhibit. 
\- It is readily seen, in the light of the familiar Poincaré theory of families of trajectories, 
that our discussion of the cases we have considered is essentially complete topologically. 
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SOME RECENT INVESTIGATIONS IN THE THEORY OF 
HYDRODYNAMIC STABILITY* 


BY 
J. R. FOOTE (Jowa State College) 
AND 
C. C. LIN (Massachusetts Institute of Technology) 


1. Introduction. In most of the existing investigations of the stability of two-dimen- 
sional parallel flows [1], two essential restrictions have been imposed: (1) at least one 
solid boundary is present, and (2) the velocity distribution is a monotone curve. It has 
been possible to treat the case of symmetrical velocity distributions through a channel 
by imposing the boundary condition at the axis of symmetry and therefore fulfilling 
the above conditions. 

There is a large class of parallel flows where no solid boundary is present. These are 
the problems of “free boundary layers’, such as the flow field in the jet, the wake, and 
the mixing region between two streams. It may be expected at once from physical 
considerations that the influence of the “viscous solutions” of the disturbance equation 
should not occur, because their influence is known to be chiefly of the nature of a 
boundary layer phenomenon in the immediate neighborhood of the solid boundary. 
In the case of the mixing region between two streams, this fact can also be easily seen 
mathematically. The exponential nature of the “viscous solutions” is well known, and 
this behavior is definitely incompatible with the condition that the disturbances must 
have finite amplitudes at large distances. The effect of viscosity must therefore be 
found from the second order approximation of the “inviscid solutions.”” The viscous 
influence is therefore restricted to the immediate neighborhood of the critical layer, 
and the general stability characteristics of such flows should therefore be very much 
different from those previously investigated. 

With these ideas in mind, Lin suggested the investigation of the stability char- 
acteristics of such flows to M. Lessen [2] (using numerical methods) and to P. Chiarulli 
[3] (using analytical methods). Lessen made detailed calculations for one profile with 
the help of Dr. L. H. Thomas of I. B. M. of New York City. The results show that the 
effect of viscosity is always to stabilize the flow even for very large Reynolds numbers, 
contrary to Heisenberg’s criterion for cases with solid boundary. The critical Reynolds 
number is found to be very low, as may be expected from the presence of the point of 
inflection in the velocity distribution and the limited extent of the effect of viscosity. 
Chiarulli obtained results for general velocity distributions of the mixing type. It appears 
from his analysis that the stabilizing effect of viscosity is probably a general character- 
istic of such flows. 

For jets and wakes it is not obvious that the above discussions apply. In the case 
of symmetrical jets and wakes, one may apply the usual type of analysis and reach 
the above conclusions concerning the origin of the viscous effects (see Sec. 2). The case 
of unsymmetrical distributions is in general more difficult, because there are usually 

*Received Aug. 21, 1949. A part of this paper is based on a thesis submitted by J. R. Foote in partial 
fulfillment of the requirements for the degree of Doctor of Philosophy at the Massachusetts Institute of 


Technology. 
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two critical points in the flow field. The purpose of the present paper is to clarify the 
nature of the asymptotic solutions in cases with two critical points. It will be shown by 
a slight extension of some work of Wasow [4] that to preserve the validity of one asymp- 
totic expression for the solutions, a path must be taken below the critical point in the 
region where the slope of the velocity curve is positive, and above the critical point in 
the region with negative slope. Application of this result to the cases of unsymmetrical 
jets and wakes again shows that the “viscous solutions” must be rejected and that the 
effect of viscosity must come from the second order approximation of the “inviscid 
solutions’. 

The remaining part of the paper is devoted to an investigation of the inviscid ap- 
proximation, with special reference to the stability of jets, wakes, and the zonal winds 
in the atmosphere. It is shown that, in all cases, a necessary and sufficient condition 
for the existence of (non-trivial) neutral and self-excited disturbances is the occurrence 
of the maximum or minimum of vorticity. This is true even in the case of zonal winds 
where the point of inflection of the profile does not have the usual significance, and the 
usual proof of sufficiency does not apply. 

It is further shown that for unsymmetrical velocity distributions the neutral dis- 
turbance may give rise to a Reynolds shear which is opposite in sign to the velocity 
gradient of the mean flow. In other words, these ordered oscillations, unlike fully de- 
veloped turbulent motions, can transmit momentum from a region of lower mean velocity 
to a region of higher mean velocity 

2. Asymptotic solutions and boundary value problems. The derivation of the linear- 
ized stability equation 

(w — ce” —ay) —w'g = - — (p'” — 2a’y’’ + a'‘y), (1) 
aR 
and the explanation of symbols can be found in [1]. We consider jet profiles extending 
to infinity in both directions. Boundary conditions for these profiles may be written 


lim ] 2 | = Fa 
yea ely) 


lim u = 0, lim v = 0. 


(2) 


ye @ y 
For the symmetric profile one may use one of these conditions together with one of the 


following: 
¢(0) _ 0, if gly) 1S odd, (2a) 


y’(0) = 0, if g(y) is even. (2b) 
It has been shown by Wasow [4] that the asymptotic solutions of-(1) may be obtained 
by formal expansions. In Sec. 6 we shall extend his result slightly and discuss how these 


solutions can be applied to profiles having more than one critical point. It will be shown 
that there are asymptotic solutions valid in a region defined as follows. The critical 


points where w(y) = c shall be denoted by y. and y , with Relw’(y.)| > O and 
Relw’(yo)] < 0. Through the point y, there are three lines C, , k = 1, 2, 3, defined by 


[i(w — c)]'” ay = 0, (3) 


ReQ.(y) = Re\ | 
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with three similar lines C¥, k = 1, 2, 3, through y (Fig. 1). Consider a path I extending 


ipn* 
i C; 





! 
/ 
C, 
I 
Gs ee. - } 
Pt i 
ye y 
Ye a 
I ies 
\ 3 
\ 
\ x 
\ a CG 
Fria. 1. 
from y = — © toy = +, and crossing only one each of the lines C, and Cf. A neigh- 


borhood of such a path T is a region of validity of the asymptotic solutions. 
The explicit expansions of the solutions are: 


g=¢ + ¢ (aR) +o)? (aR)? + -:-- 


2 = g2” + 92" (aR) + 9; (aR)? + --> 


(4) 

¢,; = exp [—(aR)'?Q.(y)]{(w — c)~”* + O(aR)””} 

¢g, = exp [+(aR)'’Q(y)]{(w — c)~* + OR) '*} 
The notation agrees with that used in [1]. The solutions g, and g, respectively become 
infinite exponentially as y ~ — © and +o. It is then clear that c; = c, = O in the 


general solution 

P = Cigi + Code + C33 + Cis 
if we apply the boundary conditions at y = +, for either the symmetric or the un- 
symmetric profile. Thus the effect of viscosity would appear to enter only through the 
higher approximations for ¢, , ¢2 . This should be discussed further for the symmetric 
profile since (2a) or (2b) may be used with one of the conditions (2). Then we are using 
only half the profile, so that yg, does not diverge and may be included: 


= CiPi + Come + Caps - 


For case (2a) the determinantal eigenvalue condition is then 


| at — , 0 
| Pir —~ APi1 $21 — Aoi 
| Fie P22 si = 0, (5) 
” ” ” | 
Pi2 P22 Ps2 


where we have put y, = —@, y2 = 0, and gf, = ¢{(y:), ¢i2 = ¢1(Y2), and so on. How- 
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ever, we shall show that the viscous effects are still not introduced by the “viscous solu- 
tions.” 
Suppose ® kg, + keg, is a solution satisfying the inviscid equation 


(w c)(y”’ — ag) — w"g = 0 (6) 


to the first approximation and satisfying the boundary condition at y, . Taking this 
linear combination of the first two columns in (5) and replacing the first column by the 
result and simplifying, we obtain the eigenvalue condition as 


Pp, = a ¥s3 (7) 
Piz 
Using (4) we find expressions for ®, and g4./¢/3 . Then (7) becomes 
1) 
P, 
@ -4 - = 
ali 
(3) 


coy? 4 Ba ne ots 
vi 7 ak 1 | uw - c)oaR [I + |. 


Neglecting terms of order 0(aR)~', (8) yields: 
o,” = 0. 


also vanishes. If we now include terms of order 


Then it follows from (6) that ®, 
Ol(akR) 3 (3) yields 


» , S 
©: ak we 
and again there is no contribution to viscosity from gy, . This process can be continued 
with the aid of the recursion formulae for g, and ¢,[ 1]. 

We can perform the same steps with case (2b) and find the condition corresponding 
to (7) as 
, 
o = of £2. 


Pa2 


Here ¢/,/{,’ has the same form as in (8), but (6) gives 


Ope w’! / 0 ) w'! “oy? 
®P, : Pp, +. a } ®P, 
Ww *} \ w—c/, 


Now w” and w c are even functions of y, so their ratio is even and the derivative 
shown is odd and hence zero at 7/2 0. Hence ®," vanishes with ®,” , and the caleu- 
lations proceed in exactly the same way. Therefore we conclude that the effect of viscosity 
enters the stability problem for the jet profile (symmetric or unsymmetric) only through the 
higher approximations to the “inviscid solutions” of the stability equation (1). This is an 
extension of the conclusion reached in the infinite jet (mixing problem) and is in marked 
contrast to other stability problems previously investigated. 

We note the following important facts: (a) the only usable solutions of (1) tend, as 
aR —o, to a fundamental system of (6), which may be obtained from (1) by formally 
passing to the limit of infinitely large Reynolds number; (b) The stability problem for 
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(1) (expressed in the eigenvalue condition) tends to the stability problem of (6) as 
ak —«. These facts give importance to a study of the inviscid equation. 

3. The inviscid equation. The inviscid equation (8) is essentially the vorticity equa- 
tion for a perfect fluid. The term —w’’(y) is known to represent the gradient of vorticity 
(1]. In other problems of stability of parallel flows (see Sec. 5), the gradient of vorticity 
also enters the problem, but it is not necessarily related to the main velocity distribution 
w(y) in this simple manner. We shall therefore consider the equation 


[w(y) — el(y” — a’y) + gly)e = 0, (9) 


where g(y) is the gradient of vorticity, and is to be replaced by —w’’(y) in the simple 
case. It is obvious that the results of Sec. 2 are applicable without modification, because 
they depend only on the coefficient w — c of ¢”’. 

Solution of (9) by means of power series gives the fundamental system* 


¢ily) = (y— y) +e? 
0 
Gye) (10) 


w'(y.) log (y — ye) 


g(y) =1tes- -— 
in the neighborhood of y = y, , and a similar system in the neighborhood of the critical 
point y = yo». The only difference lies in the path for the determination of the logarithmic 
function. The path is to be taken below the point y, , but above the point y (see Fig. 1). 
Thus, in the case of self-excited disturbances, the inviscid equation holds all along the 
real axis. But in the case of damped disturbances it fails over portions of the real axis 
in the neighborhood of y, and y, . The general solution of (9) in the neighborhood of 


y. is 
= Ag, + Be, = Ag, + Bexo —B on ¥1 log (y a Ye); (11) 


where ¢,) is the power series part of g, and A and B are complex constants. 

It is found most convenient to study the inviscid equation by examining the Reynolds 
shear stress, as was first done by Tollmien [5]. We shall now generalize his analysis from 
several points of view and then apply the results to various types of velocity distributions. 

Assuming the existence of the solutions, we shall establish three general relations. 
The actual existence problem will be discussed later. 

(i) The Reynolds shear stress is 

T= —puv = pS e'W(¢e, , ¢:), (12) 
where p is the density of the fluid, u’ and v’ are the real disturbance velocities, the averaging 
process is taken with respect to x, and W(¢, , ¢;) is the Wronskian ¢,g) — ¢gigt of the real 
and imaginary parts ¢,(y) and ¢;(y) of ¢(y) for real values of y. Tollmien discussed this 
relation only for neutral solutions with ¢; = 0. 

(ii) The rate of change of the Reynolds shear stress + along the real axis is given by 


l ne ; 
— —p= e“'cg\|eo|/|w-—c]’. (13) 
dy 2 ‘ 


*This is a slight generalization of Tollmien’s result [5). 
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It is obvious that both the above results are applicable only for parts of the real 





axis where ¢ is regular. 
(ili) In the case of the neutral disturbance, there is a discontinuity of the Reynolds shear 
at the critical point of amount 


[7] = rly. + 0) — rly. — 0) = — > mele. |'g-/| we |. (14) 


All these results can be obtained by straightforward calculations. The real disturbance 
velocities wu’ and v’ (with real values of y) are given by 
uw’ = Rey), vo’ = Re(-¥), 
with the disturbance stream function 
P(x, y, t) = oly) e* 


Thus, the average 


is easily seen to be 


- o m4 l ( me / 
Pir ~ Ori = 2; ge — gy), 
so that 
, a : : - 
HW (9, 9 ¥;) = 9 (yo — yy es (15) 


It should be noted that the conjugate signs are formal conjugates for real values of y, i.e., 
% = ¢,(y) — wily), e = gly) — ipi(y). 


From (15) one can easily derive (13) as follows. Rewriting (9) in the form 


2 gly) 
g’’ — ag g = 0, (9a) 
wy) — ec 


multiplying it by ¢ and subtracting from it the complex conjugate equation, one obtains 


ad. lo P 
(y’g —g’y) = - compces: oe 
dy * * e ¢ lw—el g 
or 
dw le |° 
Se 3 60. (16 
dy w—c| g 9) 





By combining this relation with (12), the relation (13) is obtained at once. 


1950] HYDRODYNAMIC STABILITY 271 


To obtain (14) it is only necessary to calculate the jump of W, which in turn is de- 
termined from the jump of ¢’, since ¢ is continuous across the singular point. Thus, 
from (15), 


; i.e. aca 
[W] =5 tele] —e¢le’]}, fory=y. 
From (10) and (11) it is clear that 


¢ =B, [p’) = -B- in, 
| we | 
since the angle of y — y. increases by r when y passes from y.—0 to y.+0 if w! > 0 
and decreases if w? < 0*. Thus, 


[(W] = — Tw le |. (17) 
By combining (17) and (12), (14) follows at once. 

We shall now prove the following theorem giving a necessary condition for the 
existence of self-excited and neutral disturbances. Theorems on sufficiency of conditions 
will be discussed in the next sections. 

THEOREM I. If the disturbance must satisfy homogeneous boundary conditions at 
each of two points y, and yz , the condition g(y,) = 0 for some y, in the interval (y; , ye) ts 
necessary for the existence of self-excited disturbances, and of neutral disturbances with 
w = cat some interior point of the interval (y; , Y2). 


In the case g(y) = —w'’(y), neutral disturbances with w # c in the closed interval 
(¥; , Y2) are also excluded by the condition g(y) # 0 except possibly the trivial solution with 
a = 0, ¢ = (w — ce), (which must be examined with the boundary conditions). 


PROOF: The first part can be easily proved by examining the variation of 7 for real 
values of y. At the points y, and y, , 7 vanishes. If a self-excited disturbance exists, and 
g(y) does not vanish between y, and y2 , then (13) shows that 7 varies monotonically in 
between, and hence cannot possibly fulfill the boundary conditions. 

For a neutral disturbance, 7 is a constant except for the jumps given by (14). To 
satisfy the boundary conditions, one must have 


> [7] = 0, (18) 


summation being taken for all the jumps. Since Wyin < ¢ < Wmax , there is at least one 
jump. If there is only one, we see at once that g. = 0, or alternatively ¢. = 0. We must 
show that the second alternative is impossible. Now from (10) and (11), it is clear that 
¢. = 0 implies that 


g = (w — o)f, (19) 


where f is a regular function. If we now multiply (9a) by ¢ and integrate between y, 


and y, , we obtain 


Pe Pteleltay =f gw-olsfPdy>o. (20a) 


“ws vi 


*The result, therefore also applies to the neighborhood of yo . 
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Similarly, 

| gw —c)| fl dy> 0. (20b) 
Since w — c definitely has opposite signs in the two intervals (y, , y.) and (y, , y2), it is 
obvious that we have a contradiction if the sign of g remains unchanged. 

If there is more than one point where w c, either (i) all the g,’s are zero, or (ii) 
some of the g.’s vanish or have opposite signs. The first case can be excluded by an 
argument similar to the above. The other case insures the vanishing of g(y) within 
(YW, Y 

To establish the second part of the theorem, let us note that we again have (19) 
when w ¥ cin (y¥, , y2). We may rewrite (9) as 


(w c)*f"} a’ (w ce) f + (w’ + g)(w -— of = 0. (9b 
Multiplying this equation by f and integrating between y, and y, , we obtain 
| w c)*{| f’ +a’ | f ‘t dy = | (w’’ + g)(w — cc) | f |° dy. (21 


This relation does not allow much to be inferred, except in the important case g wr’ 
In this case, we must have 


a U, f const., ne 


which proves our statement. This method of proof is well-known [6], but we include the 
details here to show why it does not work for general g(y). 

4. Wakes and jets. We shall now investigate more closely several typical problems. 
In this section, we shall study the stability of two-dimensional wakes and jets, sym- 
metrical and unsymmetrical. For these problems, g(y) —w’’(y). In the next section, 
we shall study the stability of zonal winds over a surface of revolution. 

a) Symmetrical jets. In this case, the neutral disturbance can occur only when the 
wave velocity c is equal to the flow velocity at the point of inflection, except for the 
trivial case given in Theorem I. The Reynolds shear is zero everywhere, and there is 
no transfer of energy between the main flow and the disturbance. 

For self-excited oscillations, (13) leads to a shear distribution shown in Fig. 2. One 
notes that tr > O wherever w’(y) > 0. The energy is transferred from the main flow to 
the disturbance by the Reynolds shear, as it does in fully developed turbulent motion. 
For damped disturbances, the inviscid theory is applicable only to certain parts of the 
flow field. But (13) shows clearly that there will be regions where the shear 7 has an 
opposite sign to the velocity gradient w’(y); the damped disturbance feeds energy into 
the main flow in these regions. 

We shall see that for unsymmetrical jets, 7 and w’(y) may have opposite signs even 
for the neutral disturbance. 

We have yet to prove the existence of neutral and self-excited disturbances. This 
may be done in the usual manner for the even solution satisfying the boundary condi- 
tions ¢’(0) 0 and gy > 0 as y — 0. The proof given by Friedrichs [6] can be adopted 
without change, since w(y) satisfies the boundary conditions required of g(y)..The 
existence of self-excited disturbances can be proved by following the method first used 


Ob) 


hin 


19 
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by Tollmien [5] and slightly extended by Lin [1, Part 2]. The convergence of the in- 
tegrals involved can be easily verified. 

The trivial solution ¢ = w — c for a = 0 is valid for all values of c, if one merely 
requires that the disturbance vanishes at infinity, for (u’, v’) are proportional to (ag, ¢’). 

(b) Unsymmetric jets. The existence of neutral and self-excited disturbances for such 
cases can best be seen from a continuity argument. Thus, if one decomposes the velocity 
distribution w(y) into an even part w,(y) and an odd part w,(y), one may consider it as 
a member of the family w,(y) + ew.(y). When e is zero, self-excited disturbances exist 











for 0 < @ < a, , say. When e€ changes, the eigenvalue must change in a continuous 
w 
| 
| | 
| | 
y 
| | 
| | 
| | 
| 
| 
| 
| tv | 
| | 
| 
l Cc; > 0 
| y 
| 
| 
Kia. 2. 


manner. Thus, there would still be a finite range of a for which the motion is self-excited 
and the end points of which correspond to neutral disturbances. But the range might 
shrink to zero when e is fairly large. Thus, at least for velocity distributions not far 
from the symmetrical case, we have self-excited and neutral disturbances. One notes, 
however, that ¢ = (w — c) is always a solution for a = 0. 

The shear stress distribution is very interesting in the present case. For definiteness, 
let the points where w’’ = 0 be denoted by y, and y, with y, < y, and w, > w, . For 
the neutral disturbance, it is then necessary that 


a 2 a ae 


because this is the only way to make w/’w}’ < 0. (See (14) and (18)). The distribution 
of stress 7 is as shown in Fig. 3 (dotted). The stress suddenly jumps from zero to a finite 
value at y, , stays constant, and then jumps back to zero at yo . For the self-excited dis- 
turbance, (13) gives the distribution shown by the full line of Fig. 3. The maximum 
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stress occurs at y, and the minimum at y, . One notes that in either case there is a part 
where 7 and w’(y) are of opposite signs. The net transfer of energy in the neutral case 
is, however, zero, because 


fvo 


| rw'(y) dy =r | w'(y) dy = 0. 


“We Ve 
If w, < w, , the sign of 7 is essentially negative. In either case, the shear stress in 
the case of the neutral disturbance transfers momentum from the part with lower velocity 
at the point of inflection to the part with higher velocity. 
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5. Instability of zonal winds*. We shall now apply our results to the motion of a thin 
layer of fluid over the surface of a spinning body of revolution, say, the earth. These 
results may be applied to stability of the zonal winds in the atmosphere. Discussions 
of such problems have been made by H. L. Kuo [8], largely by considering plane motions. 
With the help of present results, we can make a more thorough investigation of the 
problem. 

Let us take an orthogonal coordinate system (¢, o, z), where ¢ is the longitude, o 
is the distance from the equator along a meridian positive towards the north pole, and 
z is the distance perpendicular to the surface. Consider two-dimensional motion of a 
thin layer of fluid in the surfaces of constant z. We may put z = 0 throughout the in- 
vestigation. Let u and v be the linear velocity components relative to the earth along 


*This section was added by Lin, after the rest of the paper had been completed. For related discus- 
sions of vortex motions in the atmosphere, see a recent paper by Rossby [7]. 
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the directions of increasing ¢ and oa, and let ¢ be the component of absolute vorticity in 
the direction of the z-axis. Then the vorticity equation 


Oy , uot , of _ 
vie ee he (23) 


where r is the distance from the axis of rotation, and 


iwHabittade 
pon thie-2 cr, (24) 











> 


Fria, 4. 


with ¢, as the component of vorticity of spinning normal to the surface of the earth. 


The equation of continuity is 


ou + 2 my =O. (25) 
Oo 


0g 
Let us consider a motion slightly disturbed from a pure zonal motion 
u = Up(a), y = 0. (26) 


The perturbations wu’, v’ and ¢’ then satisfy the linearized equations 


ou’ Da 
de + ao (rv’) = 0 

(27) 
of" + Uo OF" + y! df _ 0, 


Ot r dog do 
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with 


f=i- a2. (Tuo) 
r do 


(28) 
,/ 
py lf _ 9 ay 
r \dy Oo , 
It is now possible to set up an “equivalent”? plane motion by defining the following 
quantities: 
[ do , Uo 
= y, = — 0, = 
” ‘ t io ¥* . ? 
(29) 
Poe Soe 1 eee of mt ool 
ul, = 1u’, vi = 1’, or. 
Then 
Ou% Ov% 
ss ** 
a," &, 
OS% OS% df fe 
— ee +43 Se (30) 
dt Ox, ay, 
re Ov Ou 
= = = 
sa Ox, OY, 


Thus, if we now drop the asterisks the problem discussed in Sec. 3 is reproduced with 


do 
dy (31) 


= 


The physical interpretation of the quantity 7 defined by (12) is also obvious at once. 
It is (1/27) times the total transfer of angular momentum per unit height across the 
latitude circle of radius r. 

The disturbance equation for the spherical surface has been obtained and briefly 
discussed by Kuo. 

Most of the results obtained in Secs. 3, 4 are applicable to the present case. There 
are, however, also essential differences. In the first place, the final solution has a meaning 
only when a is an integer. During the investigation of the eigenvalue problem, how- 
ever, we shall of course investigate all positive values. Another essential modification is 
the following. Since g  —w’’, the existence proof of Sec. 4 cannot be readily extended 
to the present case. It is also clear that the symmetry or asymmetry of the velocity 
distribution is not an essential feature. Rather the equality or inequality of the values 
of w at the roots of g = 0 is what distinguishes cases (a) and (b) in Sec. 4. 

The details of the investigation are also modified because the velocity distribution 
is more complicated. Typical diagrams of w(y) and g(y) are shown in Fig. 5. It is clear 
that one may consider only half of the flow (say right), and that the boundary conditions 


are 
yg = 0, or g’ = 0, at y = 0, (32) 
and 
¢g— 0, ¢g’/g = —ain the limit asy >. (33) 


| 
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But since these conditions give s = 0 at both ends, no essential modification is necessary. 
The diagrams for the distribution of 7 are similar to those shown in Sec. 4. In particular, 
one notes that the angular momentum is transferred northwards by a neutral disturbance 
if w, < w, (opposite to the case shown in Fig. 5), i.e., if the mean zonal wind at the point 
of minimum vorticity is lower than that at maximum vorticity. Otherwise, a damped 
disturbance is necessary for northward transfer of angular momentum. 


Ww 








tin 
Se — ° Yy——— P 
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The existence of self-excited and neutral disturbances will now be discussed for the 
typical velocity and vorticity distribution shown in Fig. 5. We shall consider essentially 
the neutral disturbance in the case w, = w, , because the other cases can be covered 
by arguments similar to those used in Sec. 4. For this case, with ec = w, = w, , 


=< 


K g 


is a regular non-negative function which approaches zero as y > +. Consider only 
y > 0. The equation 


ge” + re + Ke = 0, 


obtained from (9a) with —a’ replaced by \, must now be shown to have a negative 
eigenvalue. Now, the first eigenvalue is [9] 


i) 


(f? — Kf?) dy 


0 





A, = min. - 


a@ 


| f dy 
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for all functions f having proper differentiability properties, and satisfying the boundary 
condition f’(0) = 0 and making f> f° dy convergent. We want to show that for certain 
functions satisfying the boundary conditions the integral in the numerator is negative, 
for then X, is negative and our result is proved.* For this purpose, take f = e “”’. Then 


a oO 


| (f” — Kf’) dy = | (4e°y” — K)e?*™* dy. 
As e — 0, this integral approaches —J> K dy. Hence, by taking some small but finite 
value for e, the integral can be made negative, completing our proof. 

6. Validity of the asymptotic expansions. We shall now slightly extend some theorems 
of Wasow on the asymptotic solutions of the equation of hydrodynamic stability of 
two-dimensional parallel flows. Wasow established the validity of the formal expansions 
in the neighborhood of a single critical point. We shall note the modifications necessary 
to make the theorems of Wasow applicable to cases with two critical points. The pre- 
ceding sections have used the mathematical justification so derived. 

We shall retain the notation of Wasow in the main, but with x, , x* for the two 
isolated zeros; A, , A¥, k = 1, 2, 3, for the exponential-type solutions; and so on (see 


Fig. 6). 





The jet flow profiles require the assumption of two zeros for the function b (x) in 
the complex z-plane and further require that the region of validity of the asymptotic 
expressions have indefinite extent in order that the real z-axis may be used for certain 
integrals. 

The following lemma is proved by Wasow, p. 857: 





*This proof is similar to the one used by Friedrichs [6]; but the function f(y) has to be chosen differ- 


ently here. 
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LEMMA. Let T* be a subdomain of T every point of which can be reached from a fixed 
point p by a path T in T along which the rn quantities 


Re(P; — P)j), (= 1,2,---,7;7 = 1,2, ---,n), 


are monotonically non-decreasing. Then there exist r+ linearly independent solutions y = 
y-(x, A), r = (1, 2, «++ , 7) of (12) for which 


va 


y,(x, A) = ere) >> S,,(z)ra-"* + perso |, 


v=0 


where vz is an integer depending on the choice of v, , and lim », = @. 


Apparently to obtain a definite region of validity for the existence proof, Wasow 
identifies 7’ with the doubly connected domain S, and T* with S — C, . However, the 
only special property required of the domain of validity, whatever its shape may be, is 
the monotonicity condition of the Lemma. By examining the lines Re[AQ(x)] = const. 
and Re[AQ*(x)] = const. we can easily identify T and 7* with regions suitable for our 
purposes. 

In view of the equation 


Q(x) = Q(at) + Q*(a), (34) 


it is easy to show that the necessary and sufficient condition that no C, intersects any 
C*¥ is 


Re[rAQ(a*)] ¥ 0. (35) 


This condition will be assumed since it holds for any reasonable velocity profile. Then 
one sector at 2») has a common part with one sector at x¥, and by (34) the pair of curves 
bounding this domain and starting from 2¥ belong to the family 


Re{AQ(x)] = constant, (36) 


and indeed the constant is +Re[AQ(2*)]. Similarly the other pair of bounding curves 
belongs to the family 


Re{[AQ*(x)] = constant. (37) 


A set of contour lines (36) and (37) then fills the z-plane (Fig. 6) in an entirely regular 
manner except at x and 2%; i.e., through each point is one and only one curve and any 
finite region containing neither x» nor x¥ nor points of more than two sectors at x» and 
two sectors at 2* is mapped one-to-one to the corresponding region of the AQ(x) plane. 
For example, regions S, , S; , S¥, S¥ are mapped as shown in Fig. 7. In this case the 
boundaries C, , C, map into the cut C, through the origin, and C¥, C¥ map into Cf. 
Any path I of importance here passes through the common sector S; , S¥ and we may 
take it to pass also into S, and S* for definiteness. Then if we have chosen Re[AQ(x)] < 0 
in S, , on crossing C, we have Re[AQ(x)] > 0 on the remainder of I’. This holds also if 
S*¥ is used in place of S¥. It is clear from the geometrical configuration that any two 
points in the plane can be joined by a path T on which Re{AQ(z)] is monotonically non- 
decreasing, where we select the subdominant determination of the real part of AQ(zx) 
at the point p of the Lemma. This is essentially what Wasow does in choosing the point 
p, on the boundary are B, of the sector S, where Re[AQ(x)] takes on its minimum. 
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The coefficients of the basic differential equation may be taken to be such that the 
coefficients in the series (10), page 857 of [4], are regular analytic in the entire plane so 
that T can be identified with the plane, while 7* can be identified with any one of an 
infinite set of finite regions satisfying both the monotonicity property and the property 


a 
i} 


Q(x) 
AQ (xo) 





A Q(x,)=0 “i 








Fic. 7. 


that the modulus of x can be taken indefinitely large. This is sufficient for the asymptotic 
expansions since x will be fixed at some finite point while | A | ~ ©. 

Then the theorems of Wasow follow formally hy omitting the domain S from con- 
sideration and replacing it by any finite region of the above type.* In Sec. 2 we used 
this result with the following change of notation: x replaced by y; —bo(x) by w(y) — ¢; 
and {V, U, , A; , Az} by {¢1, ¢2, 3, os} respectively. 


*A different proof was given in the thesis, using only the final theorems of Wasow. 
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THE REFLECTION OF AN ELECTROMAGNETIC PLANE WAVE 
BY AN INFINITE SET OF PLATES III* 


BY 
ALBERT E. HEINS 
Carnegie Institute of Technology 


1. Introduction. In 1947, J. F. Carlson and the present writer [1], [2] investigated the 
reflection properties of an electromagnetic plane wave by an infinite set of parallel 
plates. The problem was investigated for two polarizations, that is, the cases for which 
(a) the electric vector was parallel to the edges of the parallel plates and (b) the magnetic 
vector was parallel to the edges of the plates. Let us remark in passing, that case (b) 
is equivalent to a problem in acoustics. One may see this very easily by referring to 
CH II and noting that the problem solved there is a scalar electromagnetic problem 
and one can pair electromagnetic terminology with acoustic terminology and note the 
equivalence. 











2. ke CL ot 


i 2. 








Fic. 1 


In CH I and CH II, we solved the problems subject to the fact that there be a single 
propagating mode in the ducts and a single reflected wave in free space. We shall still 
retain the first hypothesis but shall reinvestigate the second one. We shall be concerned 
with the case in which there are two reflected waves and shall find the restrictions on 
wave length, angle of stagger and angle of incidence. We shall also find the new trans- 
mission coefficient as well as the two reflection coefficients. As a check, it will be shown 
that there exists a relation between the magnitudes of these coefficients which can be 
calculated with the aid of the complex form of the Poynting vector theorem. The case 
that we consider here is the “E plane’’ case, the one which has the acoustical analogue. 

Since the formulation of this problem has been discussed in CH I and CH II, we 
only summarize the pertinent equations we used in CH II. We refer to Figure 1 for a 
view of the structure in any plane parallel to the y z plane. 

*Received Oct. 19, 1949. The material treated here is taken from an invited paper presented to the 
U.R.S.1.-I.R.E. meeting in Washington, D. C., on October 31, 1949. 
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The angle of incidence of the plane wave is 6; the angle of stagger of the parallel plates, 
which are semi-infinite in the z direction and infinite in the x direction, is a. The line 
ON 
direction of propagation of the plane wave makes with this normal. Furthermore we 
have that a r<0<aand0 < a < w/2. It is assumed here that the incident wave 


has only one component of the magnetic field, 1/,(y, z), that is, the one which is parallel 


s perpendicular to the trace of the edges and the angle # [3] is the angle which the 





to the x axis. For such an excitation, no other components of the magnetic field can be 
excited and hence the electric field may be derived from H,(y, z). If we now assume 
that the time variation of the field quantities is monochromatic, that is of the form 


exp iket), then the Maxwell equations assume the form 
| 
x E = ikH (1.1) | 
and 
V XH ikE (1.2) 
where / 2r/X and X is the free space wave-length. Here c is the velocity of light in 


units appropriate to those of A. Let us, as a matter of notation, introduce /7,(y, z) 
W(y, 2) Then (1.1) and (1.2) reduce to 


oO a : 
¥ ¥ 0 (1.3 


») 


T ; 
OY Oz 


and this is to be solved subject to the boundary condition that the z component of the 


electric field vanishes on y nd, z nd cot a. If this boundary condition is expressed 
in terms of ¥(y, z) we see that this is equivalent to the vanishing of df/dy on y nd, 
z > nd cot a, n 0, +1, --- . There are also conditions at infinity which have been 


described in CH I 
We recall from CH II, Eq. (3.2) that this boundary value problem may be formulated 
as a Wiener-Hopf integral equation for the current density on the plates. The plates 


are assumed to be perfect conductors of zero thickness. That is 


0) i in @ ¢ xp ih COS 6] 
(1.4) 


“~ t 0 i , 1 , 
} 2. | [,(z’) exp [ikup| i H, (ki (z 2’ + pg) + (y — nd pd)"}'’"| dz 


where z > 0. Here u q cos 6 + d sin @, and g d cot a, while the variable y is set 
equal to nd. Further the surface current density on the m plate is proportional to /,,(z) 


and may be expressed in terms of a single current density by the relation 
[,,(2 — mg) I,(z) exp [ikmy]. 


For all y and z we have then the following Fourier integral representation for p(y, 2). 


; wi i: 
V(y, 2 exp jsin 6 + z cos 6)] | dw J(w) exp [i(wz + kun — wgn)| 
cos {(y nd d)\(k? — w’)'’*} exp [i(ku — wg)] cos {(dn — y)(k’ — w’)'”*} 
cos [d(A w)' “| — cos (ku — wg) 
(1.5) 
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where J (w) l 


J (w) | I ,(z’) exp (—twe’) dz’. 
The path C is drawn within a certain strip of the complex w plane which we shall define 
in the following section. Here n is the largest integer in the ratio y/d. We assume here 
that k has a small positive imaginary part, and the precise meaning of small will be 
defined in Sec. 2 
2. The Fourier transform of the kernel of equation (1.4). The Fourier Transform 
of the kernel of the integral equation (1.4) is 


"th wy" pM exp lip(ku — wg) + | p | d{k? — w’}'”). (2.1) 
We choose here that branch of (k? — w*)'” which reduces to k when w = 0. The series 
in (2.1) is not in a particularly useful form since it does not exhibit the singularities of 
the function it represents, save for the possible branch points at w = +k. We say possible 
branch points advisedly, since in closing this sum, we shall find that they disappear 
and the only singularities possessed by (2.1) are two infinite sequences of poles. In any 


event, each term of the series has branch points at w = -+-k, so that each is regular in 
the strip Ymk < Ymw < Ymk. We observe further that the series (2.1) diverges when 

ku — wa + dik’? — w}' ? = 2nr (2.2a) 
or 

ku — wg — d{k? — w*}'” = 2ne (2.2b) 
n 0, +1, 42, --- . From Eqs. (2.2a) and (2.2b) it is clear that the singularities of 
(2.1) are intimately related to the values of a, 0, d and k which we use. For n = 0, we 
find that w may be k cos (2a — 6) = o, or k cos 6 = o, . Since a — r < 6 < a and 


0 <a < w/2, it is clear that Ymo,. < Ymo, . From this we see that the series converges 
at best in a strip Yme, < Ymw < Ymo, and this is so, provided that there are no other 
singularities given by Eqs. (2.2a) and (2.2b) which are in a smaller strip. 

A next possible root which may be in the strip {mo, < Smw < Yme, comes from 
qs. (2.2a) and (2.2b) by putting n = +1. Let us assume that the angle 8 = 2/2 + 


6 — a satisfies the inequality 0 < 8B < 2/2. Upon solving Eqs. (2.2a) and (2.2b) for 
n l we get 

w, = (2r/d) sin af — cos a(1l — ku/2mr) + Al (2.3a) 
and 

w, = (2r/d) sin a[— cos a(l — ku/2xr) — Aj (2.3b) 


where A is given by the expression 
A = [(dk/2x)? — sin? a(1 — ky/2m)*)]'” 
and the positive root is understood for A. The roots w, and w, are written in this form 


under the assumption that A is real for k real. Observe that w, satisfies Eq. (2.2a) while 
w, satisfies Eq. (2.2b) for n = 1. 
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We now write A, A for k real and assume further that A does not vanish. Then 


we may expand w, and w, for QYink e sufficiently small (where now e is small enough 


to insure the convergence of the expansion) and we obtain 

Ynw, (e sin a)/d][u cosa + (d°k, + 2rp sin’ a ku? sin’ a) /2rA, | 
and 

Ymnw, l(e sin @) d\{u COS a (d°k, + 2rpu sin’ a key sin’ a) 2rA,| 


if we only retain first powers of e. In terms of the angle 8, these expressions may be 


simplified to read 


Ye e(cos asin B t sin a(dk, COS” B + 2r sin @ sin B) 2rA, | 
and 
Ymnw e[cos a sin B sin a(dk, cos B + 2m sin asin B)/27A,]. 


Clearly, under the restrictions we have imposed on @ and B, Ymw, < Ymw, so that it 


remains to compare the relative orders of magnitude of Yinw, , Yinw, , Yio, and Ymo;, . 


We shall now show that Siw, < Qme, . This implies that 
cos a sin B sin a(dk, cos’ B + 27 sin asin B)/2rA, < cos 6 
or 
cos asin B sin a(dk, cos’ B 2r sin asin B)/27rA, < sin (a + £) 


From this, we get immediately 


sin a(dk, cos’ B + 227 sin asin B)/27A, < sina cos B. 


Since the left side of this inequality is always negative for a and @ positive and acute, 
while the right side is positive, it is clear that Qnw, < Qino, . We have similarly that 
Yinor < Yinw, , Yo, < Yimw, , and Yinw, < Yino, . If we now combine all of the 


inequalit ies we get 


Yinw, < Yo. < Qinw < Ymo, < Ym, 


Hence the series (2.1) converges:in the strip 

Yintos << Ynw < Ymo, 
provided now that there are no other roots of Eqs. (2.2a) or (2.2b) in this strip. Such 
will be the case, if we limit dk, properly. Let us suppose for the present that such is the 


case 
One observes at this point, that if A, is real for real k, that is for k = k, , then the 


following inequality is satisfied by /, , d, @ and a. We have here 


dk, sin a 


2r 1 + cos (0 a) 


But in order that one mode propagate in the ducts, 


0 < (dk,)/2x < 3 
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Hence, subject to the condition that 0 < 8 < 2/2, the strip of regularity of (2.1) is as 
we have given it, provided 
sin a dk, 1 


< 


( - - ° 
jak OX ES te 2 


This inequality on (dk,)/2m assures us of the existence of at least the two propagating 
reflected waves in free space. Since we desire only two such waves, we should put down 
conditions which insure that the remaining A, forn = —1, +2, +3, --- be imaginary 
for k, real. In the first place, the remaining roots of Eqs. (2.2a) and (2.2b) are 
|(2an sin a)/d|[— cosa(1 — ky/2rn) + iA,] 
where 
A, = sin’ a(1 — ku/2xn)? — (dk)’/4n’x’]'”’. 
A,, will vanish if 
(dk)/2ne = — sina/[{1 — cos (@ — a)] n= +1, +2,--:- (2.4) 
or 
(dk)/2ne = sin a/|1 + cos (6 — a)] n= +1,42,:--. (2.5) 
But for n negative only, the roots 
dk/2nxe = — sina/{1 — cos (@ — a)] n= —l, —2,-:- (2.6) 
satisfy the inequality dk/2x > 0, while for n positive, only the roots 
dk/2nxr = sin a/{1 + cos (@ — a)] n=1,2,--- (2.7) 


satisfy the same inequality. 

From the two sequence of roots (2.6) and (2.7) we can now form the condition that 
there be only two reflected waves. One reflected wave arises from the root w = a, . 
This we have seen in CH II corresponds to the reflected wave which makes an angle 
with respect to ON equal numerically to that of the incident wave. The other one corre- 
sponds to w w, . In order that a third reflected wave not be present 


(dk,)/2r < (2 sin a)/[1 + cos (@ — a)] or sin a/|1 — cos (@ — a)]. 
If dk, satisfies this inequality, it satisfies the inequality 


0 <a- 6 <f 2 
(dk,)/2r < (nsin a)/[1 + cos (6 — a)]} 
n = 3, 4, 


as well as 
O<a-—-6<7/2 


(dk,)/2e < (nsin a)/[1 — cos (@ — a)] 
n = 2,3, - 


We have then the following condition for two reflected waves. 


0 < (sin a)/[1 + cos (@ — a)] < (dk,)/2r < {} or (2sin a)/[1 + cos (@ — a)] 


or sin a/[1 — cos (@ — a)]}} 
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The term on the right side is to be interpreted as the smaller of the three terms. For 
example if sin 6 = .25, this inequality reads 


0 < (4sin a)/5 < (dk,)/2x < {4 or (4sin a)/3} 
On the other hand, if sin 8 = .5, this inequality reads 
0 < (2sin a)/3 < (dk,)/2r < {4 or (4sin a)/3}. 


It is clear that we can find conditions for three or more reflected waves by arguments 
similar to those that we have just used, but we shall not pursue this detail any further. 
One final remark is in order. Suppose that +/2 < a —@ < 7m. In this case the inequality 


reads 


0 < (sin a)/[1 — cos (6 — a)] < (dk,)/2x < {4 or (2sina)/[1 — cos (6 — a)] 


or sin a/[1 + cos (@ — a)]} 
with the same interpretation lent to the right side of the inequality as we did for the 
case 0 < a — 6 < w/2. 


3. The solution of equation (1.4). We define Eq. (1.4) for all z to read 


F(z) = Y(z) 
(3.1) 


“ ; Z. | T,(z’) exp [kup] “ Hs" [k{(z — 2’ + pg)? + (y — nd — d)*}"”"] dz’ 
eed ae 07 


where F(z) vanishes identically for positive z. Further ¥,(z) = zk sin @ exp [zkz cos 6] for 
positive z and vanishes identically for negative z. In view of the restrictions we have 


placed on dk, , 6 and a, we have that 
F(z) = Ofexp (to.z)| + Olexp (7w.z)] (3.2) 


for z— — ~. In order to see this, we rewrite Eq. (3.1) in Fourier integral representation as 


F(z) = (z) — 1 | dw J(w) exp [i(wz + kun — wgn)] 
dor 


“¢ 


{k? — w’}'”* sin [d{k? — w*}"”7] 


" cos [d{k? — w’}'”] — cos (ku — wg) — 


Since J,(z) has the same z dependence as the magnetic field, and since only the lowest 
mode propagates in the ducts, J(w) is regular in the lower half plane §mw < Ym. 
Further since there are two reflected waves, the path C is drawn in the strip {mo, < 
Y¥mw < 3Ymo,{4]. In order to find F(z) for z — — ©, we need only calculate the residues 
of the poles in the lower half plane {mw < Simo, . Upon doing so, we find that the domi- 
nant terms are of the form given in Eq. (3.2); all other terms which have been omitted 
are small compared to these two. This incidently tells us that the unilateral Fourier 
transform of F(z) is regular in the upper half plane {mw > Sime, , provided we assume 
the integrability of F(z) for z negative and finite. This may be verified with the solution 
of the integral equation (3.1). The Fourier transform of Y(z) is regular in the lower 
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half plane {mw < Ymo, . Finally the transform of the kernel (2.1) may be written in 
closed form as 


{k? = w?}? sin [d{k? _ w'}'/?] _ K_(w) 


cos [d{k® — w*}'”] — cos [ku — wg] K.(w) 





and this is regular in the strip Y¥mo, < Ymw < Ymo, . There is then a common strip of 
analyticity to all transforms involved and it is thus permissable to apply the Fourier 
transform to Eq. (3.1) to get 





_ksin @  K_(w)J(w) | 


D(w) wo — @, 2K .(w) 


If we separate terms as we did in CH II we get 


2k sin 0K .(¢,) 
J(u) = ooo 
(w — o,)K_(w) 
and 
ksin 0K .(w) — K, 
Diy = £9 GE Meo — Efe 
(w — o,)K,(w) 
The integral function of separation vanishes identically in this case. The terms K_(w) 
and K.(w) are now given by the following expressions. 


K_(w) = N_(w)/D_(w) and K.(w) = D.(w)/N.(w) 


where 
N_(w) fexp [r(w)}{k — w} I] [{1 — (kd/nm)?}'? + iwd/nx] exp [—iwd/nz] 
n=} 
N.(w) = {2d(k + w) 
x exp [—x(w)]/(d? + g’)} T] ({1 — (kd/nm)*}'? — iwd/nx] exp [iwd/nz] 
n=1 
fr ) 
D_(w) = (w — ,)§ [] [4. — ®,] exp [(ku — wg + iwd)/2nr + i(x/2 — a)] 

Xx { [] (4, + *,] exp [(ku — wg — iwd)/2nx — i(x/2 — 2} 

xX {(%, — A) exp [(ku — wg + iwd)/2x — i(r/2 — a)]} 
and 


( 1 
D.(w) = (w - a8 I] [A, + 2¥,] exp [(ku — wg — twd)/2nx — i(nr/2 — «)} 


( « 
*< I] [A, — 7W,] exp [(ku — wg + twd)/2nr + i(x/2 — 2} 


x {(v, + A) exp [(ku — wg — twd)/2r + i(x/2 — a)]}. 
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Here 
A, fsin* a(1 — ky/2nmr)? — (kd)?/4n?9?}'”” n= —1, +2,--- 
A {(kd)*/4mr° — sin’ a(1 — ky/2m)*}’ 
Y, [cos a(l — ku/2nm)|] + (wd ese a)/2nr n = +1, +2, 
and 
x (w) [twd/2r|[(a — 2/2) cota — In (2sina)]. 
The subscripts — and + on A, N or D indicate that the function is regular [5] in the 
lower or upper half plane, as the case may be. For example D_(w) is regular in the lower 
half plane SYinw < Ymo, , N_(w) is regular in the lower half plane ¥mw < QYim, ete. 
I ‘ 
4. The investigation of the far fields. We have found J(w), the Fourier transform of 
[,(z). If we now use Eq. (1.5), we have an integral representation for ¥(y, z). From this 


— = 


representation we can determine y(y, z) for z © (in the ducts) and y(y, z) for z > 
(in free space). For z > 0, we close the path C in the upper plane [6] and obtain 


isin 0K .(o;) 


v(y, 2) oT KR .(b) [1 — exp (ku — tkg)] exp [i(kz + kun — kgn)] 
sank — O))INGK 


+ Of(cos ry/d) exp [—z{x/d° — key] 


Hence for z ~@, ¥(y, z) is asymptotic to 


isin 6K .(0;) = en : 
Dd(k — 0) KB [1 — exp (tku — ikg)] exp [i(kz + kun — kgn)] 


and from this we can find the transmission coefficient in the duct. It is 


isin 6K .(o;) sid 
7 = —sin K .(o, : (1 — exp (iku — tkg)]. 


2d(k — o,)K.(k) 
For k real, | 7 |° reduces to 


_— (% — 4) (¥: + a) sin’ (a — 6) 
Vv, + A/,\v, — A/,, [cos 0/2][sin’ (a — 0/2)] 


where the subscripts k and o, on the first two fractions indicate that w has been evaluated 
at these values. 
For z < 0, we close the path C in the lower half plane and obtain 
¥(y, 2) = exp [ik(y sin @ + z cos 6)] 
4 J(o2) exp [ik{y sin (2a — 6) +2 cos (2a — §)}] 
2d[cot (2a — 0) — cot a] 


J(w2) {hE — wi}! exp [ifwee + (KY = wi)" 1) 
4r(cesc a)A 
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terms which attenuate exponentially for z + — ©. The term (k* — w;)'” may be written 
as [2x/d] [A cos a — sin® a(1 — ku/2r)]. Hence the first reflection coefficient is 

R, = J(oe2)/2d{cot (2a — @) — cot al 
while the second reflection coefficient is 

R, = —J(w.)(k? — w>)'”/4rA ese a. 
For k real 


= sin’ 0 K.(o,) |? 
R. = > eee 7 


(o2 — o,) sin’ (2a — @) | K’(a>) 





where the prime in K.,(w) denotes differentiation with respect to the argument w. Upon 
simplifying further, we get finally for k real 


R, = [cot’ (a — 6/2)]|t in’ 6 2i| ¥ =4| [z +4) | 


Furthermore for k real 


R. |? = 4k*[sin® a|[tan” 6/2]|sin* (a =O) +4) k oo | 
—_ [w. — o,|[w. — 02] Vv, — Aj..LAE — wey 


2 


The values which we have given for | 7 |, | R, |” and | R, |? have been derived subject 
to the conditions 





( . 
(: we sin a <$ dk t < 1 = _2 sin a sie sin a \ 
- 1 + sin B lo" 1+sinB. 1—sing 


(b) 0 < B < 2/2. 


The modification for —7/2 < 8 < 0 is quite direct. We observe that for this case, 
r/2 <a— 6 < x. The requirement for two reflected waves for this interval of a — @ is 





(oe 0O< __sin @ kd < {dor — 2 sin a Sina _\ 
7 1 —sinB 2r 2 _ sin B” a + sin B)° 
In view of the requirement on a — @, it is clear that the inequalities (a) and (c) are 


equivalent. There are a number of modifications in the expressions for | R, |*, | R2 |’, and 
T | for the range —2/2 < 8 < 0, but we shall not pursue the matter further. 
5. Conservation of power flow. As a check for our calculations of | T | ,| R, | and 
R, | we shall now see what information we can derive from the Poynting vector theorem. 
In the first place, the power flow from free space into any duct is (per unit length in the 
x direction) 


a / (2 a3 x H*)-n dy 


where n is the normal vector into the parallel plate region. In terms of y, this may be 
written as 


oe ~d 
Re -~< me | | cos a 4 — sina at lye dy. (4.1) 


“0 
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Here y = z tan a, that is we are integrating over the opening of a parallel plate region. 
Hence upon simplifying the parenthesis and integrating over y, we get 

d esta . : ‘ | ” 

> {kK sin (a — 0) —k | Fe, [sin (a — 6) — 2rA | R, \"/d]. 

This represents the average flow of power into a given duct per unit length in the x 
direction. This in turn is to be balanced by the average flow of power far down the duct, 
that isasz—-o, nd < y < (n + 1)d. The magnetic field in the duct for z ~~, nd < 
y « n+ 1)d is 


¥(y,z) = T', exp [zkz]. 


Hence the power flow down any duct for z © is | 7' |*d/2 (recall that | 7’, | is inde- 
pendent of n, the duct number). Hence conservation of power flow demands that 


oh ae i esc a{k sin (a — 0) —k| R, |’ sin (a — 6) — 297A | R, |*/d] (4.2) 


Upon substituting in | R, |° and | R, |° we find that Eq. (4.2) is an identity, thus verifying 
the conservation theorem. 

6. Some concluding remarks. If we try to find conditions for higher reflected waves, 
that is the cases corresponding to n 1, —2, --- orn = 2,3, --- , we find that the 
conditions for a given number of reflected waves become fairly complex. On the other 
hand, a study of the zero of the A,’s and their ordering will show us how each reflected 
wave enters. These results will depend on @ and 8. The case 8B = 0, has been excluded 
from our present discussion, but it offers no difficulties. In this case, we get higher re- 
flected waves in pairs. The result is not surprising since the incident wave normal and 
the normal to the trace of the edges of the plates are parallel. 

l’or the case of a single reflected wave, we noted that the reflected wave has a normal 
whose direction is independent of *. In the case of two reflected waves, the first reflected 
wave has a normal whose direction is independent of k, but the second reflected wave has 
a direction which is dependent on k. We shall present a discussion of the numerical 
results elsewhere. 

Finally, we observe that | R, |°, | R.|° and | 7'|° have not been normalized. If we 
refer back to Eq. (4.2) and divide through by & ese a sin (a — 6), the normalized first 
and second reflection coefficients, as well as transmission coefficients for the case dis- 


cussed, are respectively 
RI, 27rA ese (a 6) | R, |°/d and T |’ sin a esc (a — 6). 


For the normalized coefficients, all magnitudes are numerically less than unity. 
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4. At this point we realize that to get more than two reflected waves, the inequalities in sec. 2 have to 
be modified. As we carry on this modification we find that an indefinite number of reflected waves can- 
not enter. For example, if a = x/2, then the inequality for two reflected waves reads 


1 dk . fi 2 1 
0< hj < eb <befe 


for 0 < 6 < x/2. This is impossible since sin @ < 1. Hence there is only one reflected wave in this case. 
The formulation we gave in CH I makes no assumptions as to the form of ¥(y, z) to the left of the parallel 
plates. The convergence study in sec 2, of the present paper gives us conditions for one, two, etc., re- 
flected waves. 

5. See CH I, p. 321 for some remarks regarding the present use of the term regular. 

6. See CH I, sec. 4. 
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—NOTES— 


A VARIATIONAL PROBLEM FOR THE ROOTS 
OF A CUBIC EQUATION* 


(A contribution to the theory of servomechanisms) 
By HANS BUCKNER (Minden, Germany) 


1. The problem. Let ¢ denote the time, u(t) a function of time and v(¢) the displace- 
ment of a servomechanism, which is to simulate u(t). As regards the error 


n(t) = v(t) — u(t) 
many servomechanisms obey the equation 
7” + a,n”’ i a7’ a a3n = —x!"" — au” (1) 


with constant coefficients a, , --- a, . The quality of the servomechanism depends on 
the roots of the algebraic equation 


2 +a,2° + az +a; = 0. (2) 
These may be denoted by 2, , x2 , x; . We consider the quantity 
E = Min (—R2, ’ —Rzr2 ’ —Rz;), 


where $tz, means the real part of 2; . Then — > O means stability, and the greater , the 
faster the servomechanism follows sudden changes of u(t). A positive — appears if and 
only if the well known Hurwitz-conditions are s:.tisfied. These are 

a, > 0, a, > 0, a,a, > a; > 0. 


If a, and a, are regarded as fixed, and a; as variable, ¢ is a function of a; . An interesting 
problem both from a mathematical and a technical point of view, is to find Max &. There 
is a simple formula for it, which does not seem to be as well known as its merits. 


2. The solution. Introducing the classic notations 


a; a> 
3 ’ q=Q-— 





p=a@- 


D = —4p* — 279’, 


we may represent the roots by 


os a; _ q iI is 1/2 a os ¢ _ 1 me ai 
wh 7} 6s 9 a se 1/2 a £ q_ - = wn” 

t= bt a 9 + jg (-3D) | + 4 — 5 — jg (8D) ; (3) 
ee a, om q 1 7 si 1/3 E q ™ 2 - ai 

r= m4 | 9 + 7g (-3D) | + 9 ~ 1g (~3P) 


*Received March 28, 1949. 
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where 
! 1 
») 
The quantity p depends on a, and a, only. 
, t are continuous functions of a, 
a,a,. In the first case one root vanishes, 


The roots x, , # Hence the same holds for €. 


The quantity & vanishes for a, 0 and for a, 
‘A aa 
t-i(a,) “~ represents two roots with vanishing real part 


in the latter a i 
For 0 < a, a,d, we have — > O. Therefore & as a function of a, has a positive maxi 


mum in the interval (0, a,a,). 


The criterion for multiple roots is D = 0. For D > 0 three different real roots exist; 


for D) < 0 one pair of conjugate complex roots exists. The roots 2, and x, of the system 


(3) give these complex roots, while x, is real. 


Phe only case where 2 may vanish is for p 0. Vanishing takes place for two values 


of a, which we denote by (a B) and (a + £8). We can write 
a0 ‘\1 
(1, Qs 2a; ' p 
a ee 3 2 = 
" 9 9 
” Poe | a 
All roots are continuously differentiable functions of ay in all intervals which do not 
contain (a 8) or (a + B). Starting from the well-known formulae 
Yr, +a +a Gs « 
ri t+ ata ay 2a 
UV, tt +a a, + 3a,a, 1 
we find 
dv, l 
da, (vy Ke)(U, “i 
dr. I 
; (4) 
da, (al Uy) Xs) 
Ax. l 
da. (7 L)(L5 Le) 


Now, let. a, be a value for which D > 0. Then all roots are real and different from one 


another. Hence —£ coincides with one root, say 2; , for a sufficient small neighbourhood 
of a, . From (4) it follows that 
di dx, 


' ~ ), 
da, day 


in this neighbourhood. Thus, a maximum of & can exist within the interval (0, a,a2) 
only if D < 0. 

Next, we consider a value a, and a neighbourhood of it in which D < 0. In this 
case we may write 


o.. o}.. ay l q I ‘ 1/2 l 9" 1/2 
Mar, Ma 3 Rt ig | 3D) Ot eee ee 18 | 3D) 
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and, furthermore 
§ = Min (—z, , —2z): 
The criterion for 
z, = Rez, R 


is q 0. This means that a, = a. 
If it is assumed that a, does not coincide with a, — is represented either by —2z3 or 
by —Mz, for a sufficient small neighbourhood of a, . From (4) it follows that 
lt dx, 
= —KR—! <0 
da, da, 
with a suitable subscript 7. This holds for the neighbourhood. Therefore a maximum 
does not exist if D < O and q# 0. 
Summing up, we may state that the maximum of ¢ is given either by D = 0 or by 
q 0,D < 0. 
The following three cases may occur. 
1) p> 0: We find D < 0. Max € is determined by gq = 0. We have 
a, . A, 7 


as = 


Max é = 3 ' 9 97 a: « 


2) p 0: We find D 0,q = 0 for a, = @ and 


9 
a ’ A» rr 
Max é = 3 for a; =~ — oa - 
‘ e | 


3) p <0: It is obvious that D > 0 for g = 0. The maximum of & is obtained for 
D = 0. The roots then are given by 


1/3 


1/3 
_. 2) 
Zs 3 - 2( 4 , 


For a, a + B, we have q = +8 and 
3 1/3 
t= : — (2) for a, =a-+ B, 


a B 1/3 
_ ; = 2(8) for a, =a — B. 


-~ 


‘Tr 


We see that a, = a + 8 gives the greater value for ¢. As the maximum of € exists, at 

least one of the abscissae a, = a + 8 must be in the interval (0, a,a.). Since D is negative 

for dy = a@,a, (where complex roots exist), the abscissa a; = a + surely is in the interval. 
Therefore, we may state that 

2 


9 (3)'"(—p)*”. 


Q,Q, _ 2 
3 27 


1 l /2 ‘ 
Max ¢ : — 3)” (—p)'” for a; = a; + 
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Obviously we can combine all cases and write 


f 1, 
Maxé = (4 - - 73 (—p)' 23 
AX § 3 (3) 72 (—p) 


This is obtained for 


PERIODIC MOTIONS OF A NON-LINEAR DYNAMIC SYSTEM* 
By H. SERBIN (Purdue University) 


1. Introduction. It is a matter of common knowledge that many problems of dy- 
namics are non-linear in character and that the linearization of such problems is an 
expediency adopted in the face of mathematical difficulties attending the non-linearity. 
Fortunately, the linearized treatment gives useful results in a number of problems. On 
the other hand, there are cases where the essential characteristics of the phenomenon 
are altered by the assumption of linearity. 

In the absence of general methods for handling non-linear problems, it is desirable 
to investigate typical non-linear phenomena to as great detail as possible with the 
intention of arriving at general methods of analysis. 

One such problem is van der Pol’s equation 


- il. a eee £2) 
x —el—27)rx +2 =0, (: =U’ = ae (1.1) 


This defines the one-dimensional motion of a particle of unit mass suspended on a spring 
of unit stiffness and subjected to negative damping at low amplitudes and positive 
damping at high amplitudes. 

The existence of periodic solutions of an equation of more general form 


a + f(xjx + g(x) = 0 (1.2) 


was considered by Liénard, Ref. 1, who proved, for g(x) = x and under fairly general 
conditions on f(x) (see Ref. 2), that there exists a periodic solution. The proof used in 
Ref. 1 was reproduced in Ref. 3 for g(x) ¥ z. 

In the present paper, the problem is considered from a different point of view. The 
periodicity of the solution is shown to exist for a wider class of differential equations 
(1.2) which can best be described and compared as in the following table: 


az “Zz 


F(a) = | f(x) dx G(x) = | g(x) dx f(x), g(x) 


Lefschetz F(2) =@ G(o“) > 0 Symmetric 
Present paper F(o) > 0 G(o) =o No symmetry required 
F(o©) =o G(o) > 0 No symmetry required 


*Received July 6, 1949. 
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In addition, the case in which F() and G(@) are positive and both finite is dis- 
cussed. It is shown that an example can be constructed in which a periodic solution 
does not exist. 

A number of inequalities relating to the initial energy and amplitude of motion are 
derived on the basis of which one can, in some cases, calculate upper and lower bounds 
of the amplitude of the periodic motion without solving the differential equation. This 
should be useful in application to practical problems. It is conceivable, in some instances, 
that an average of the upper and lower bounds will define the amplitude with sufficient 
accuracy for the purpose required. 

Since the problem is non-linear, one must expect such inequalities to be non-linear, 
for instance, algebraic. The solution of these inequalities is carried out conveniently by 
graphical methods. This must not be compared with graphical methods of solution of 
the differential equation. It is conceivable that for a given degree of accuracy, examples 
may be constructed of the type (1.2) in which the graphical method of solution by 
isoclines would give results incorrect both qualitatively and quantitatively, showing a 
periodic solution where none exists. 

It should be pointed out that the assumption of symmetry in f(z) and g(x) made 
heretofore is not consistent with the nature of the problem. If these functions are de- 
termined empirically, then one will not have symmetry in the mathematical sense if 
x and —2 represent two different positions of the system. One is accustomed, in small 
deflections, to replacing an empirical curve which is nearly parabolic at the origin by an 
exact parabola. However, for large deflections, this is not permissible. Any lack of sym- 
metry vitiates the validity of the theorem of periodicity as known heretofore. 

A mechanical system illustrating the non-symmetric case is afforded by a pendulum 
in which friction acts proportional to velocity. If one adds a device which acts through 
a limited range of motion of the pendulum at the bottom of the swing in such a way 
as to sense the velocity and then to force the motion proportional to the velocity over 
the friction, then one would have the features characteristic of Eq. (1.2). If the constant 
of proportionality is different for z and —2z, the asymmetric case will exist. 

2. Qualitative character of the solution. Consider the differential equation (1.2) and 
define 


F(a) = [ f(a) ae, 


G(x) = [ g(x) dx, (2a) 


H = 3x? + G(2). 
Make the assumptions: 


(A) f(x) and g(x) are continuous for « 2 0, 
(B) f(x) < OforO0 <xr<24,, 

(C) f(x) > Oforz>x2z,, 

(D) F(o) > 0, 

(E) g(x) > O for x > 0, 

(F) F(@e)G(oe) =o, 


Then F(x) has only one positive zero, x = 22. 








298 NOTES [Vol. VIII, No. 3 


Since 


dH 


= x(x + g(x)], 


it follows from Eq. (1.2) that 


dH 


= —f(x)x°, (2.1) 


which shows that H increases with time for x < x, and decreases for x > 2, . Dividing 
both sides of Eq. (2.1) by x and simplifying, 


Of = ¥ fla)(2H — @))”, wa 


where the sign is minus or plus according as the velocity is plus or minus. 

Suppose H plotted versus x in an (x, H)-plane (Fig. 1). It is seen that H = G, the 
equality implying that dH/dx = 0. At a point where H > G, the right side of Eq. (2.2) 
satisfies the Lipschitz condition in H. Since f(x) and G(x) are continuous, it follows 
from the fundamental existence theorem on the solution of ordinary differential equations 
that there is one and only one solution passing through (z, H) for H > G. Referring to 


Fig. 1, there are four distinct parts of the trajectory: 


Part x f = 
j ym 

01 + _ + 
._ + fm 
a6 - + + 
56 _ — _ 


The subscripts 0, 1, etc., will indicate the value of a variable at the corresponding 
point of the trajectory. 

It is clear that H, , H. , --- , He , 23; are monotone, continuous, functions of Hp . 
The proof given in Ref. 3 shows, in addition, that H, — Hp» is a monotone increasing 
function of zx, . 

The analysis in the following pages is directed toward showing that H, — H, , re- 
garded as a function of H, on the interval 0 Ss H, <~@, changes sign. Because of the 
monotone character of H, — H, , it will follow that there is one and only one trajectory 
for which H, = Hi. 

3. Behaviour of the solution for 0 < x < xz, . Consider only those trajectories for 
which 7, = zx, . Define 


L(G) = | (2H — @)'? dH, 


where the subscript 7j indicates integration along the trajectory in the sense 7 — j. If 
G(x) < G’(x) on the part 7j of the trajectory, then J;;(G) < I,;(@’). Since 
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1;; = ¥ [ f(x) dx = +F(F; — F;,), 
it follows that’ : 
I,,(0) = (2H,)'? — (2H,)'? < 11,(@) 


(3.1) 
= —F, < 1,,(G,) = [2(H, — G,)]'” — [2(4, — G,)]'”, 
I,,(G,) = [2(H, — G,)]'” — [2(H. — G,)]'” < 1(@) 
(3.2) 
= —F, < 1,,(G,) = [2(H, — G,)]'” — [2(H, — G,)]'”. 
From (3.1), it follows that 
{(2H,]'? + F,}? < 2H, < {[2(H, — G,)]'” + F,}? + 24,, 
and from (3.2) that 
([2(H, — G,))'” + F,}? > 2G, < 2H, < ([2(H, a G,)]'” i F,}? + 2G, ° 
The latter two inequalities yield 
0 < 2H, — H.) < —2F,{[2H,]'” — [2(H, — G,)]'”} 


a —4FG, —4F,G, 
— 12H]? + (21H, — G»]'”} * GH)” * 

The smallest value of x; for which the derivation of (3.3) is applicable is 7; = 22 ; 
for this choice of x, , it follows from the first part of (3.3) that H, = G, > H, . Hence 
when H, exceeds G, , x; will be larger than z, and (3.3) is valid. 

A similar argument applied to the trajectory 456 yields 


0 < 2(H, — H,) < —2F,{[2H.]"? — [2(H; — G,)]"”} 


eG AG, — 
{(2H,]'”” + [2(H, — G,)]'”} {(2H,)'”* + F,}’ 
where the last inequality in (3.4) is based on the assumptions H,; — G, > O and 


([2H,]'* + F, > 0. However, when H, = H, , these are both satisfied if H, is chosen so 
large that 
(2H,)'” = (2G,)'” — F,. (3.5) 


Fx Tr, 


(2H;)'” = (2H,)'” — Iss(0) > (2H,)'” + F, = (2H_)'”? + F, = (2G,)'” > 0. 


Therefore, assuming H, = H, and (3.5), one has from Eqs. (3.3) and (3.4) 
2F,G, 2F,.G, 


H, — Hy < (Hy — He) + He - HY) < — Gay — GH, 
_— oo : 
= — 27.) oH) + ((2H,)” rs =) 


1The following analysis, which leads to inequality (3.3), sharper than the author’s original result, 


is due to the reviewer. 





(3.6) 
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4. Behaviour of the solution for x = x, . Since J,.(@) = J,;(@) = F; and (H — G)s2 > 
(H — Gu, 


H,-H,=| aH>| aH=H,-M, 


“32 “43 
(4.1) 
H, — H, = (H, — H;) + (AH; — H,) < 2(H2 — Hs). 
Also, 
I5.(G2) = [2(H, — G.)]'” — [2(H; — G,)]'? < Is2(G) 
(4.2) 
= F, < I3.(G;) = [2(H. — G;)]'” < [2(H. — H,)]'”. 
From the first part of (4.2), it follows that 
2(H, — G.) < [2(H; — G.) + F3J’, 
(4.3) 
2(H, — H;) < F; + 2F,([2(H; — G,.)]'” = F; + 2F;[2(G; — G,.)]’””. 
Combining with (4.1), we obtain 
H, — H, < F; + 2F,[2(G, — G,)]'”. (4.4) 
Using the second part of (4.2), we have 
LF? < H, — H, < Fi + 2F,[2(G, — G,)]'”. (4.5) 
Using (3.6) and (4.5), we can now state: 
Theorem 1. When conditions (A)-(E) are satisfied and H, is so large that 
(2H,)'” = (2G.)'” — F,, (4.6) 
f 4 1 } 
7m > —4F,.G) 5-43 + cme) s 
Py 2 —4F, “any” ((2H,)'” + F,])’ sade 


then t, > x, and H, < Hy. 

As H, increases, H, increases. If H, — H, >, then, from (4.5), F(~)-G(e) =a. 
On the other hand, F(@) =© implies H, — H, —~&, by (4.5). If G(@) =, then 
H.(> G,) approaches infinity as 7, ©. Hence by (3.3), H, ©. Therefore condition 
(F) is equivalent to: 


H,-@ as 37 @. 


Now suppose that the differential equation (1.2) is defined for x < 0 and that condi- 
tions corresponding to (A)-(F) are satisfied for x < 0; specifically, assume that: 


(A’) f(x), g(x) are continuous for z S 0, 
(B’) f(z) < 0 for Ge <2 <0, 
(C’) f(x) > 0 for s<s, 

(D’) —F(—-) > 0, 

(E’) —g(x) > 0 for 2 <= 0G, 

(F’) —F(—~)G(—o) =o, 
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If one replaces x by —&, eq. (1.2) becomes 


é + (OE + g€) = 0, 


where 
g.i(f) = —g(-8), 
Then 


nt E 
10,90 FhO=[ hed Geoe=[ aa 


satisfy conditions (A)-(F), although with x, replaced by —2{ . 

Now consider in the region x = 0 a trajectory starting with initial energy H, , re- 
turning to « = 0 with energy H, , continuing into the region x < 0 (with initial energy 
H,), and returning to x = O with energy H,,. . Because of assumption (F), H, and H,. 
are monotone increasing functions of H, in the interval O S Hy, < Hj. . According to 
Theorem 1, H, < H, for sufficiently large H, . If Hs were unbounded as H, —~, then, 
applying Theorem 1 to the region = — x = 0, H, may be chosen so large that H,. < He. 
Hence for large Hy , Hi. < Hy. 

Now consider motions of small amplitude. When z; S x. , Ho < H,.A similar 
result holds for negative x. Therefore for small Hy , H,. > Ho. 

Therefore ,, is a continuous, monotone-increasing function of H, for0 S H, <@ 


such that H,. — Hp» changes sign. For some H, , Hi. = Ho , showing that there is a 
periodic solution of eq. (1.2). It is shown in Ref. 3 that H, — H, is monotone decreasing. 
Then H,. — H, and H, — H, are also monotone-decreasing and the periodic solution 


is unique. 

Theorem 2. Jf Eq. (1.2) satisfies conditions (A)-(F) and (A’)-(F’), there is a unique 
periodic solution. 

5. A counter-example. Let us denote a trajectory recurrent when H, = H, . It is 
shown in Ref. 3 that if f(z) is even and g(x) odd relative to x = 0, a periodic char- 
acteristic is recurrent. It will now be shown that when condition (F) is not satisfied, 
there is, in general, no recurrent trajectory and therefore no periodic solution. 

Consider the family of equations (1.2) defined by 


IIA 


x: 


2) 


F(z) = F(a), 0 


IIA 


x 


>xr, O<eS1 


| 
~ 
x 
= 
& 

V 


For every «, 


G, = H, S G(o) <om. (5.1) 

Now for every choice of H, satisfying (5.1), consider the trajectory 456 (Fig. 1) 
and that trajectory 012 for which H, = H, of the afore-mentioned part 456. Then 
H, — H, is a positive-valued, continuous function of H, on the closed interval (5.1); 


H, — H, therefore attains a minimum A > 0 for some Hi, . 
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If there were a recurrent solution for (1.2) for every e,0 < ¢ S 1, then from (4.4), 


O0<A Ss H, as H, < [F.(xs)]? + 2F (x3) {2[G(x3) aaa G,]}'”" 


e[F(@)]° + 2eF(@){2[G(@) — G,]}'”*. 


A 


By choosing ¢ small enough, the inequality may be violated. Hence for such an e, 
H,~ H,. 
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Fia. 1. Qualitative characteristics of the trajectory. Fig. 2. Graphical solution of upper bounds. 


6. Numerical bounds. It is useful, for application to specific problems, to develop 
numerical bounds to the amplitude x, of recurrent motion. Since x; > 2x, , a lower bound 
is x, . The following set of inequalities is useful, in some cases, for obtaining an upper 
bound to xz; and H 

For x3; > 22, one has, using (3.3) and (4.3) 


— ‘ 2F',G, 
G, < H, + (MH. — Ho) < Hy - (2H,)' 5 (6.1) 
H, < (H, — H;) + H; < £ + F,[2(G, — G,)]'? + G, . (6.2) 


A trajectory satisfying inequalities (4.6), (4.7), (6.1), and (6.2) has the property H, > H, 
and its maximum amplitude x, and initial ordinate H, furnish upper bounds to the 
corresponding characteristics of the periodic motion. 

Example. Take f(x) = —3(1 — 2°), g(x) = z. 
Then, 


F(z) = 1(_, oo ; ), G(z) = . ' z, = }, 2, = 3", 


F, = G, _ 


The graphical solution of the inequalities is shown in Fig. 2. The curve of zx; vs. Hy 
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lies between the curves labelled (6.1) and (6.2); the intersection with curve (4.7) 
defines a point whose coordinates are suitable as upper bounds. Referring to Fig. 2, 
it is seen that x, at point A is an upper bound to the amplitude and H, at point B is 
an upper bound to the initial energy. Thus 


1.73 < 23 < 2.50 


and the average 2.11 of the bounds is in error by less than 23%. 
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QUENCHING STRESSES IN TRANSPARENT ISOTROPIC MEDIA 
AND THE PHOTOELASTIC METHOD* 


By R. C. O'ROURKE anp A. W. SAENZ (University of Michigan) 


Introduction. It is well known that when a transparent non-crystalline solid, such 
as glass, is heated to a uniform temperature 7’, and then rapidly quenched in a bath 
at temperature 7)(7, > 7) there results a non-uniform stress distribution. Depending 
on 7’, , one can divide these stresses into two distinct classes. If 7’, is sufficiently lower 
than the softening temperature of the material, then the stresses are transient, but if 
T, is sufficiently high (say 550° C for lime glass) then the quenching stresses are not 
transient but, on the contrary, remain permanently set into the glass. The latter are 
referred to in the literature as quenching or residual stresses. As is well known, the 
existence of such residual stresses is made manifest when the object is examined in 
polarized light, and by violent explosive characteristics of quenched objects when cut. f 
Each distribution of such stresses is characterized by a definite double refraction pat- 
tern. In 1841, F. E. Neumann’ developed a general mathematical theory of the double 
refraction of light in non-uniformly heated isotropic solids. In turn the problem was 
studied theoretically by such men as Maxwell,’ and Lord Rayleigh.* 

The purpose of this paper is three-fold: 

1) To develop a mathematical theory of residual stresses based on a simple model 


*Received July 11, 1949. 

tFor details on the background of the subject see E. G. Coker and L. N. G. Filon, Treatise on photo- 
elasticity, Cambridge, 1931, §§ 332 and 333. 

1F. E. Neumann, Die Gesetze der Doppelbrechnung des Lichtes in comprimierten oder ungleichférmig 
unkrystallinischen Kérpern, Abh. K6énigl. Acad. Wiss. Berlin, Part II, pp. 1-254 (1841). 

2J. C. Maxwell, On the equilibrium of elastic solids, Roy. Soc. Edin. Trans. 20, 87-120 (1853). 

3Lord Rayleigh, On the stresses in solid bodies due to unequal heating and on the double refraction 
resulting therefrom, VI-1, pp. 169-178, 1901; see also Arch. Neerl. (II) 5, 32-42, (1900) and Collected 


Papers, vol. 4, p. 502 
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of the quenching process. This theory is essentially a careful reinterpretation of the 
classical theory of thermal stresses, and is implicit in the work of Lord Rayleigh.* 

(2) To justify, at least in a preliminary way, the use of the photoelastic method in 
the study of quenching stresses with particular emphasis on glass. 

(3) To recognize the problem of determining the distribution of stresses photo- 
elastically in symmetrically quenched spheres and cylinders as a simple problem in 
integral equations. 

1. Photoelastic preliminaries. In this section the necessary photoelastic formulas 
used in the study of residual stresses in spheres and cylinders will be developed. The 
Maxwell-Neumann stress-optic law* states that the secondary principal stresses in the 
plane of the wave front are related to the retardation dR, suffered in traveling a distance 


dy, as follows: 
dR = C(P — Q) dy (1.01) 


provided the direction of the principal stresses in the plane of the wavefront does not 
vary along the light path (Fig. 1). It is desired to calculate the integrated retardation 
suffered by a plane monochromatic wave in passing through symmetrically quenched 
spheres and cylinders. The wavefronts are parallel to the z-z plane (Fig. 1). For the 
integrated retardation one can write: 


R(n) = c| (P — Q) dy (1.02) 


(C is the stress-optic coefficient in Brewsters; y is measured in mm; (P — Q) is the differ- 
ence of the two principal stresses in the plane of the wavefront). 

In the case of a long cylinder, whose surface is kept at 7’, , the stress tensor in cy- 
lindrical coordinates assumes the following diagonal form: 


g,(p) 0 0 
|| 
g;; || = 0 o4(p) 0 |i, (1.03) 
0 0 a (p) || 


where o, , o, , and a, are functions of p only. In order to obtain the principal stresses 
in the plane of the wavefront, one must transform the stress tensor (1.03) to appropriate 
cartesian axes (See Fig. 1) by means of a rotation about the z axis through an angle ¢, 
given by 
cos sin ¢ 0 || 
R,@¢) = || —sng cos d 0 | (1.04) 
0 0 1 || 


which transforms the stress tensor to: 


‘For the basic photoelastic concepts, see Coker and Filon, loc. cit., Chapter III, especially paragraph 
3 and 10. 
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1950] 
cos’ ¢o, + sin’ do, sin ¢ cos ¢(¢, — a4) 0 
Oz Try Trz 
| Tye C, Ty, || = | sin ¢ cos ¢(¢, — G4) sin’ do, + cos’ day 0 (1.05) 
| 
Tez Tay o, | | | 
0 0 Oo. | 
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Fic. 1. Normal incidence on cylinder. 


From (1.05) one readily identifies the principal stresses in the plane z — z of the wave- 
front as 


P=go@,, 
(1.06) 


Q = cos’ do, + sin’ de, . 
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Reversing roles of P and Q merely changes the sign of the retardation. This is equivalent 
to inverting the Babinet Compensator.” Thus the total retardation becomes: 
a | 


R(n) = C | {o. — cos go, — sin’ do,} dy. (1.07) 


> 


Now, this integral can be simplified because of the inherent symmetry involved in the 
study of spheres and cylinders. One expresses cos’ ¢, sin” g, and y in terms of £ (Fig. 1) 
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Fig. 2. Normal incidence on sphere. 


and p, and notes the fact that o, , ¢, and o, are functions of p only. Then equation (1.07) 


becomes: 


f | pdp 


R@) = 2C | {( — 6) + (6, — ¢,)° (1.08) 


2 2 2\172 
pi\p —€&) 
This equation is of fundamental importance for the study of residual stresses in 
cylinders and spheres. It can be studied in two different ways, both of which are treated 
in detail below. The function #*(£), (an asterisk will signify that the function is obtained 


’Coker and Filon, loc. cit. paragraphs 1 and 36. 
I gray 
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experimentally), is easily obtained with a suitable Babinet compensator. (1.08) reduces 
to an Abel integral equation with the aid of appropriate substitutions, thereby enabling 
one to give an exact solution for each of the three stresses, ¢, , 7, , ¢. . (See Section ITI). 
The second way of treating (1.08) is to employ a model of the quenching process (Sect. 
II), and calculate o, , o, , and o, , thereby finding a theoretical retardation pattern 
R(é) which can then be compared with experiment. The case of the sphere lends itself 
to an analogous treatment. One measures R*(£) along the equatorial plane (See Fig. 2) 
where no rotation in the plane of the wavefront occurs. The stress tensor is diagonal in 
spherical coordinates r, 3, g, for a sphere quenched in the same way as the cylinder 
above. It assumes the form: 


| o,(r) 0 0 | 
o:; || =|] 0 o9(r) Ol; a=%. (1.09) 
0 0 o,(r) 


The retardation integral for the sphere becomes: 





. ab Ir 
RG) = 208 | (e, — 0) Tan * (1.10) 


2. Applicability of photoelastic method in the study of quenching stresses. In Sec. 3 
a method will be developed which reduces the calculation of residual stresses in sym- 
metrically quenched spheres and cylinders to the inversion of Abel’s integral equation. 
For the cylinder, this method will depend upon the assumption that residual stresses 
are elastic stresses; this assumption is implied in the “sum rule.” The “sum rule” is 
the relation o,(p) = o,(p) + o,4(p) which is characteristic of thermoelastic stresses in 
symmetrically quenched circular cylinders (hollow and solid). The necessity for making 
this assumption lies in the fact that the method being proposed uses relative retardations 
R*(£) (Babinet Compensator) instead of absolute retardations (Mach-Zehnder inter- 
ferometer). If absolute retardations are measured, then one can invert the appropriate 
Abel integral equations of Sec. 3 without the ‘sum rule” ansatz. (This will be reported 
on elsewhere.) However, for practical purposes, it is much more expedient to measure 
relative retardations since, as will be shown in Sec. 3, one obtains the complete distri- 
bution of residual stresses in symmetrically quenched cylinders from a single photo- 
graph of R*(é). 

It is therefore necessary to give some justification for the “‘sum rule” ansatz in the 
case of residual stresses of the type being considered. There are at least four exact 
methods for showing the elastic nature of these residual stresses, two of which will be 
explained now. First of all, Eq. (1.08) can be written in the following form: 

_p dp 


ob 
R*(é) = 2C | OOF — By» (2.01) 


since equilibrium demands that 
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Now if one uses the “sum rule” and the equilibrium equation 


+ te) = 0 
dp PT 5; = 


Eq. (2.01) becomes 


rae | 
R*(~) = 2C | — (p'a,) an ° (2.02) 


From (2.02), using Dirichlet’s principle for inverting orders of integration and the 
fact that o,(a) = o,(b) = 0, one easily derives the interesting relation 


| R*(&) dé = 0 (2.03) 


For those residual stress patterns which have been studied, (2.03) is very accurately 
satisfied. Secondly, it is of further interest to see if one can construct a model of the 
quenching stresses from which a theoretical expression for R(€) can be computed and 
compared with the experimental pattern R*(é). This is indeed possible and will be briefly 
discussed here. First of all, one assumes the existence of a “‘state of ease’’ at a temperature 
T*, such that for 7 > T* all stresses in the body are negligibly small compared to the 
stresses for T <« T* (of order 10° mm’? for lime glass). Secondly, one neglects the 
temperature variation of the elastic and thermal constants. The striking agreement 
obtained so far between theory and experiment seems to indicate that the model is 
essentially correct. Many experimental investigations have revealed the fact that the 
properties of amorphous materials, such as glass, change rather abruptly at a certain 
critical temperature.” The proposed model merely idealizes this fact by making the 
elastic constants change discontinuously at the critical temperature 7*. 

These assumptions readily lead to the following rule: 


To calculate the residual stresses in a body (A) quenched from the uniform 
temperature 7', to the temperature 7’, of a bath (7', > 7), one calculates 
the thermal stresses at time t = ¢* (7 = 7%) in an elastic body (B), with 
the same elastic and thermal constants as (A), whose initial temperature is 
T, and whose surface beginning at t = 0 is maintained at 7, . The stresses 
in (B) at ¢ = ¢* are the negative of the residual stresses at t © in (A). 2.04) 


Since the problems of thermoelastic stresses in spheres and cylinders have been dealt 
with thoroughly in the literature, it is not necessary to discuss them here. Essentially, 
one takes the pertinent solutions (say from Timoshenko’), replaces ¢ by ¢*, changes the 
signs of the stresses, according to (2.04), and substitutes into (1.08) and (1.10) to obtain 
the theoretical retardation patterns. The theoretical retardation patterns for cylinders 
and spheres have been calculated in this fashion and are given below: 

‘R. Houwink, Elasticity, plasticity, and structure of matter, Cambridge, 1940, Chap. 3 and 6. 

6S. Timoshenko, Theory of elasticity, New York, 1934, Chap. 11. 
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Solid Cylinder (radius b): 





\ _ _ 4CaK(T, — To), 5 exp (—xéat*) 
ade toe 2 eT EH) Py 
(2.05) 


—)**'(21 + 1): -(&,b)***?# 
fo foo |= 4*** bk + 214+ 21)! 





Me: 
Ms 
att 


R(n; ) = —J,(&,.b)(1 — d*)'? + 


where 





JorsiQgsb) {_1- 2k + 0x ! 
d lQk+D+1°2k+)-1 


are positive roots of J,(t) = O and A = &/b (Fig. 1) 

In practice ¢* is sufficiently large so that only the “‘shape factor” R(1, \) need be 
calculated (see Fig. 3): 

(B) Hollow Cylinder: (a < p < b) 


“(J = rer, 


g,, 

















RQ) = —2C =—— >> exp (—xét*) + Aa {20 Z(G.) — 2abZ, (6.0) 
‘ b* — a’ 
(—)*****(21 + 1)-(E,b)***? Zorar(Aknb) 
2.06 
7 e > > 4**'k\(2l + k + 1)! r 7 
, - nart-e) 1 —Ad° 2k + Ir? 
alti. tare + eS} 
where 
_r ee —— JSo(Ea) wz / Ps 7 
A, = 5 (T; T En Te. b) — Ja) {bZ,(E,b) aZ ,(E,a) } 


Z (x) = Jo(x)NolE.b) — No(x) Job) 


£, are roots of Z,(ta) = 0 


(C) Sphere (radius b) 


a YA 3\ 1/2 os 
RQ) = - sah. (T, - 7»)(©) -b- } (—)'n exp (— ney ea a) TF d), (2.07) 
n=1 
where | 
*(2n)*-'?k! 





R(n; ) = >! “Ok +)! (1 — n7)***? J, 45/2(ne) 

The leading terms R(1; \) for (2.05) and (2.07) are plotted in figures 3 and 4 re- 
spectively. The agreement between these “first approximations” R(1; A) and the corre- 
sponding experimental patterns R*(\) obtained so far furnishes strong argument for the 
applicability of the photoelastic method in the determination of the type of quenching 
stresses under study here. 
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3. Photoelastic determination of quenching stresses in spheres and cylinders. 
Equation (2.02) is Abel’s integral equation so that the problem of finding the distribution 


of residual stresses in symmetrically quenched cylinders reduces to the inversion of this 
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integral equation, which is well known. The complete solution for cylinders can be 


written in the form: 


. 4 ER*(E) dé 
o,\ p) 3 | ~3 5173 » 
TU p (e — p) 
d ; ‘ 
o4(p) (pa,(p)), (3.01) 
, dp 
a(p) = o,(p) + o4(p). 


For the case of spheres, one immediately recognizes (1.10) as Abel’s integral equation 
without the necessity of a “sum rule” (which does not exist for the corresponding 


thermoelastic stresses). A similar inversion yields: 


o,(r) = ’ 


f° — r*) 


2 [ R*(&) dé 
TU E( 
(3.02) 
r a 


g(r) = o,(7) = 5 ad TC, 
«4 ar 





mn 
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Thus, the entire problem reduces to quadratures. R*(£) is an analytic function of 
the light-path variable n (or (b” — £*)'”*). Thus one can write 
R*(t) = do a,(v? — &)'” (3.03) 
k~1 
If (3.03) is substituted into (3.01) the quadratures can be performed and one finds 
for cylinders: 
’ (1/2) < r'(k + 1) 2 2\ (k+1)/2 
7o,(p) = — = a (BF — ptt? 3. 
p*a,(0) Ff 2 Reree ~" (3.04) 


Correspondingly for spheres equation (2.07) suggests the following ansatz: 


ce) 





R*(é€) = > cb" = gf)?’ (3.05) 
and (3.02) yields 
I i(k + 1) 2 2) (k+1)/2 ¢ 
ar) = — (1/2 ee —r 3.06 
i 7610/2) Dares gy U - 1) esas 


The coefficients a, , c, are easily obtained by fitting R*(é) by the method of least 
squares. 

From the above results, it is seen that the method in this section rests on the ‘“‘sum 
rule’ for the case of concentric and solid cylinders, while in that of spheres no such 
assumption has to be made. It is thus clear that equations (3.01) and (3.02) are inde- 
pendent of the detailed nature of the boundary conditions of the problem, and that they 
are therefore applicable when large radiation losses are present which cannot be accounted 
for by the simple Newton law of cooling, a situation which would preclude the use of the 
method of Sec. 2, from a practical standpoint. 

The results obtained so far employing the procedures of Sec. 3 are of technical im- 
portance in locating exposed tension regions (i.e., unstable regions) in quenched cylinders 
and spheres. In such work one desires to know the sign of the boundary stresses, their 
absolute magnitude being of secondary importance. One could also gain rough ideas 
concerning residual stresses in non-transparent polycrystalline materials, such as steel, 
by simply changing the elastic and thermal constants appropriately and using ideas of 


similarity. 
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STEFAN-LIKE PROBLEMS*' 
By G. W. EVANS II, E. ISAACSON, ann J. K. L. MacDONALD? (New York University) 


Introduction. Consider the recrystallization of an infinite metal slab, as for example, 
the change from a to 6 crystals when iron is heated through the temperature 1643°F. 
Accompanying this change of crystalline form is a latent heat of recrystallization. If 
the slab is originally of uniform crystalline structure and is heated from the front face 
through the critical temperature, new crystals appear at the heated face. If the front 
face is planar and is heated uniformly and the back face is insulated, the interface 
between the two types of crystals is planar. For the determination of stresses occurring 
in the metal, it would be necessary to determine the history of the interface as it moves 
inward. To determine the history (a space-time curve) of this interface, it is necessary 
to solve the heat equation, in this case, a linear homogeneous parabolic equation, with 
appropriate boundary conditions and with a suitable discontinuity in the first order 
space derivatives across the unknown interface. 

In each of the several problems to be discussed here, we will consider the recrystalli- 


zation of a metal slab of uniform thickness, L, bounded by the planes x = 0 and x = L. 
We assume that a uniform heat source covers the front face, x = 0, and that a perfect 
insulator covers the back face, x = L, if L is not infinite. Thus, the problems will involve 


only one space dimension. Furthermore, the coefficient of thermal diffusivity, a’, is 
assumed to be a constant for a given type of crystal, and recrystallization is assumed 
to take place at a temperature u = 0. Since in each problem the object is to find the 
propagation of the interface through the slab, the choice of the critical temperature as 
zero in no way affects the results. 

To determine the type of discontinuity across the interface, assume that the recrys- 
tallization has already occurred to a distance x where 0 < x < L. The subscript 1 
after a symbol will be used to indicate that the symbol refers to the region of the original 
crystalline form and the subscript 2 to indicate that the symbol refers to the recrystallized 
region. The equation of heat balance across the interface as it moves a distance dx then 


is, after division by dt, 


op, H2x(t) = (kote, — kyu,,) |, = x(b, (1) 
where x = 2(t) is the space-time curve along which the interface moves, x (t) = dx/dt, 
p is the density, H is the latent heat of recrystallization, and the k; are the coefficients 
of thermal conductivity. The notation used is du;/dx = u,;, fort = 1, 2. The temperature 


distribution in each region is given by the differential equation 


Uz, = aui,, 


*Received Oct. 14, 1949. 

1a) The term “Problem of Stefan” has recently been used for problems similar to the ones treated 
here by L. I. Rubinstein, A. B. Datsev, and others. 

b) The authors express their indebtedness to the members of the Institute of Mathematics and 
Mechanics and especially to Professors R. Courant, B. Friedman, and J. Keller. 

2Deceased. The practical results of this paper were obtained prior to the untimely death of our inspir- 
ing co-author J. K.L. MacDonald. The surviving authors shall endeavor to complete, in the mathematical 


sense, the material covered in this note. 
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1. The first problem to be considered will illustrate a method of solution for all. 
In this problem the slab has any thickness L, 0 < L < @. The slab is considered to have 
been heated uniformly to the critical temperature, 0, and then a constant heat source 
is applied to the front face of the slab. This may be represented by the following mathe- 


matical problem: 


Find x = x(t) where the temperature distribution satisfies the equations® 
Ul, = A Us: for0 < x < 2(d), (3) 
u=0 for allz > 2x(d) (4) 
with the boundary conditions 
A x(t) = u,[x(0), ¢] (5) 
x(0) = 0 (6) 
u(O, t) = g, (7) 


where g is a constant. The condition (5) is a consequence of equation (1) where A = pH/k. 
We obtain a power series representation for « = x(t) about t = 0 by assuming u 
expressible in a power series about x = 0, t = 0 for all u in the region 0 < zx < z(t), Le., 
any discontinuities in the derivatives of u are assumed to occur only on crossing the 
curve x = x(t). 
By differentiation of (4) and (5) along the curve x = x(t) and of (3) in the region 
0 < x < x(t) and use of (7), all derivatives of z(t) may be determined at the origin. 


a(t) may then be written as the following Taylor series: 


, g I ; q° 2 5 g° 3 
nO 94 - gee THR 
(8) 
Sa «, Ww » 
oe ee Oe et = vee 
4! QA’ 5! a®A® 


The coefficients of the series have not been calculated beyond those given. They may 
be found from the following relations: 


Letu= >> a,,a't' and z= >) C,t', then, 


i,7=0 i=0 


> (i+ Naser i( > Bh ct’) t=A DL (nt+ Cat’, (9) 


t,7=( n=O 


) a.i( > ct’) i = 0, (10) 


i,7=0 n=0 


1 (Wgt) (11) 


ant > -s ° co a: a; ,; 
ao (i+ 1)\i+2) °°" 
with a,,, = g and C, = 0. 
’We have adopted the convention of using the letter u without subscripts when the temperature is 
identically 0 in the region 1 (i.e., to the right of the discontinuity). 
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In the event that the heat input, g, is no longer a constant but an analytic function 


of time, g = g(t), we find in a similar way the interface curve to be: 


( Cs 3 ( 4 Cs 5 
HQ Ct tte te +e mee, (12) 
a o vo: 
where 
(0) 
¢, = 
A 
— gq’ (0) q (0) 
~ @ A® A” 
— 59°) _ 6 g (0 g°(0) 1 g_ (0) 
aA’ a A A 
= 519°(0) _ 75 g (0) q'(0) 4 10 g (0)9°O) 
a’ A’ ao A a A 
15 [a (0 }* g(0) g (0) 
* a A” _ 
rp — 827 9°) _ 1438 g (0)9"(0) 4 195 g (0)g*(0) 
a A a’ A’ a’ A® 
525 [g (0)]’9°(0) 15g (0) gO) 


15 ig (0 60 g (0) g (0) g(O) g_ (0) 
a eee eee a. 


a’ A a’ A 


The series of Eq. (12) reduces to that of Eq. (8) for g(t) a constant, g(0) = g. 
2. To generalize the treatment of Sec. 1 to cover the case of non-analytic boundary 
and initial conditions, we employ Laplace transformations. Furthermore, we believe that 


the questions of uniqueness and existence will be answered by these methods. In this 


part, we derive the integral equation which is satisfied by the interface curve x = x(t). 
To this end, the discontinuity curve is considered in the form ¢t = f(x), and for greater 
generality we consider the heat input, g = g(t), to be a function of time. The problem 

is to solve the differential equation 
u,(z, t) = a’u,,(z, t) for t > f(x) (13) 

with the boundary conditions 

u(z, t) = 0, fort < f(z) (14) 
u,(0, t) = g(t), (15) 
f0) = 0, (16) 
. A fa 
u,[z, f(z)] = =, (17) 





mn 
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where f’ = df/dz. 
The Laplace transform for u(z, t) is defined as 
uz, 3) = [ u(x, te“! dt = [ u(x, t)e*' dt. 
/0 “ f(z) 


Application of the transform to (13) furnishes the transformed differential equation 


E - 5 thal adil (18) 


with the boundary conditions: 
uF(0, s) = gs), Lim {u [(z, s)} = 0. 


The solution of Eq. (18) is 


pd 1/2 1/2 
We.) = = | exp {—sf(t) | {- a la—é i + exp f. os (x + o} dé 
“@\5) “0 a 
ag™(s) " \ 
— I eee ui, 
s a 


and the inverse transform of this equation is 


oT | a oo 5a 
ute, 0 = > | ine — FO” [ox fF t— f®) 


+ en aaa fe ao 


_, [fivt=-9 a} 
aj (nn) exp {=e dr. 


Once t = f(x) is known, Eq. (19) gives the temperature distribution in that region which 
has been recrystallized. Now, applying the condition (14) in the form u[z, f(x)] = 0, 
Eq. (19) becomes 


. , a ‘a sittin i ote lo ben et 


f —(x + &)° \) 
+ exp 4 E ; 
Pe 4a? [fC t)— f()] 7 £(#) 


Equation (20) is an integral equation for the determination of t = f(x). To solve Eq. 
(20), we return to the designation x = z(t). To do this let t = f(x), 7 = f(é), x = x(0), 
and £ = x(r). This change is permissible if t = f(x) is monotonically increasing. After 
making these substitutions in the right hand side of Eq. (20), replace t — + by 7’ and 
replace 7’ by 7 in the right hand side. Equation (20) becomes 


2 [ a ath exp {= iC) Db a = A [ st 9 | exp = tO nt = si 


[x(t) + a(t — DP i a. (21) 





(20) 








4a’ r 


b- 
+ exp § 


7T>0 
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A power series solution for x = x(t) of equation (21) may be determined as follows: 
Let 

a(t) = Cit+ Ct? +---, (22) 
and assume g(t) a constant, g. The coefficients C, and C, are determined by using (22) 
to expand the integrands of equation (21) in powers of ¢ (i.e., 2(¢ — 7) is expanded about 
7 = 0). We find that in the limit, as ¢ tends to zero, equation (20) becomes 


r= 





Now, replacing x(t) by gt/A + C.t° and z(t — r) by 


/ oz ( 
a(it— r) = 2) — re ()O+r ia 


in equation (21) and expanding about ¢t = 0, Eq. (21) becomes 
1 8 , 1 





9,2,-1/2 
4gt'’? — = + g 3 = 4gt'” 
, aA A’a ; 
16 A , Qq°n'”” ll gays 
=~ 0,0? —- 145 , 
3 % aA + B 364 zl 
which yields 
1 g° | 
C= — —s 
‘ 2a°A* 


It is noticed that the coefficients determined in this manner are the same as the first 
two of equation (8). Similarly, when g(t) is an analytic function of t, the first two co- 
efficients, C, and C, , of equation (22) have been determined and agree with those of 
equation (12). The determination of the remaining coefficients by this method becomes 
extremely tedious. 

3. Another problem which involves the determination of the curve of recrystallization 


would be to prescribe an initial temperature distribution in a metal slab of infinite 
thickness (L =o). The problem is then to determine the curve x = x(t) where the 
temperature distribution is given by 
Us, (x, t) = agle,,(z, 2) for0 < x < x(d), (23) 
u[a(t), t] = 0, (24) 
u;,(z, t) = aju,,,(z, b) for z(t) <x <o@ (25) 


with the boundary conditions 


kotto,{x(t), t] — kyw,,[x(t), t] = pH x(t) = Bx (0), (26) 
u(x, 0) = o(z) where d(x) < O and ¢(0) = 0, (27) 
(28) 


U2(0, t) = g(t). 


,) 


it 
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In this problem g(t) is assumed to be the same as in the previous problem and ¢(z) is 
assumed to be expressible in a Taylor expansion about x = 0 with an infinite radius of 
convergence. The method of solution is again the same as that of Part I and the series 
representing x(t) is given below through the second power of t. 


I , 
s(f)) = B |k.g(0) — k,o’(0)]t 


+ - {h| - oF (k2g(0) — kye'(0))” 
. (29) 


+ v0) | m | £0 (kag(0) — k,6’(0)) 


This problem is reducible to the second problem of Part I by taking ¢(7) = 0, and the 
coefficients in Eq. (29) then agree with those of Eq. (12). 

Comments. The boundary condition (4), « = 0 for x > x(t), used in the first two 
problems is not as specialized as may first appear. In most practical cases the metal is 
initially at a uniform temperature, say room temperature. If a metal slab, initially at 
a uniform temperature, is heated at one face and is insulated at the other, the temperature 
gradient remains zero for sufficiently small L. That is, the temperature rises uniformly 
until the critical temperature is reached providing the slab is thin. This result is seen 
if one lets u = U > 0 be the critical temperature and considers the slab initially to have 
the temperature distribution u(z, 0) = 0. One can then show that the temperature 
distribution is independent of the space coordinate as long as u(x, t) < U. The tempera- 
ture distribution in the slab of thickness L is given by the differential equation 


Ur = Urs (30) 
with the boundary conditions 
u.(0, 2) = gd), (31) 
u(L, t) = 0, (32) 
u(x, 0) = 0. (33) 


The Laplace transform of (30) is the equation 
nm s o 
u(x, 8s) = —; u(z, 8), 
Qa 


the solution of which is 


1 O(s) cosh (L — x)s'’?/ 
u(x, 8) = oe owls . ie Poe s (34) 





for the transposed boundary conditions (31) and (32). Since the slab is considered thin, 
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1/2 . 1/2 . . . ° * 
replace cosh (L — x) s ‘“/a and sinh Ls’’"/a by the first terms of their expansions giving 


? a” Or 
u-\x,s) = — jj g-(s). (35) 
$8ly 


In (35), w(x, s) is the transform of u(z, t), g2(s) of g(t) and —a’*/(sL) is the transform 
of —a’/L; therefore, u(x, t) may be written as the following convolution 
ia a” 
uz, = - / TL IY) dy. 

This expression shows u(x, t) to be a function only of the variable ¢, and is therefore a 
constant with respect to 2. 

Discussion of the curves of equations (8), (12), and (29). a) Figure 1 shows the degree 
of agreement of the first five sums of the series of equation (8) where x,(t) represents 








t 
0-54 
X, 
x 
2 X3,Xs5 
0-44 
ra 
> x, 
2 o34 
< 
Ww 
= 0-24 
= 
0-14 
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O-1 0-2 o-3 * 
THICKNESS IN FEET 


Fia. 1. 


the first term on the right-hand side of (8), x.(t) represents the first plus the second, 
z;(t) represents the first three terms, etc. The calculations for this figure are based on 
the following data: k = 25 BTU/(hr.)(sq. ft.)(°F/ft.), c, = 14 BTU/(lb.)(°F), p = 
480 lbs./cu. ft., H = —6.76 BTU/Ib., and g = —65 BTU/(sq. ft.)(hr.). This data is 
approximately that for a cast iron. 

b) Figure 2 shows a comparison of the results of Eqs. (8), (12), and (29). For this com- 


parison the boundary and initial conditions are taken in the following way. The curve 


for Eq. (29) was drawn using the following data: p = 480 Ibs./cu. ft., H = —6.76 
BTU/lb., k, = 25 BTU/(hr.)(sq. ft.)(°F/ft.), k, = 26 BTU./(hr.)(sq. ft.)(°F/ft.), 
¢,, = .12 BTU/(lb.)(°F), c,, = .14 BTU/(lb.)(°F). g(t) was chosen to be —117 e7' 


BTU/(sq. ft.)(hr.) and ¢(x) to be —50 tan™'z °F. This choice of g(t) is such that the 
heat added per unit time decreases with increasing time, while ¢(z) has the property 


that the initial temperature at x = 0 is the critical temperature and the temperature 
decreases with increasing x until it reaches 25r°F below critical at x = ~. For the curve 


of Eq. (12) the data is acquired from that for Eq. (29) in the same manner that the 


solution (12) may be obtained from the solution (29), i.e., kz = k, ¢, = Cp, o(z) = 








r 
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¢(0) = 0. The data for the curve of equation (8) is deduced from that for Eq. (12) 
similarly, i.e., g(t) = g(0) = —117 BTU/(sq. ft.)(hr.). Since only two terms of the 
solution (29) have been determined, all curves are drawn using only the first two terms 
of each series. The straight line shown in Figure 2 represents the time it would take a 
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THICKNESS IN INCHES 


Fic. 2. 


slab of thickness x at the critical temperature to recrystallize uniformly with a constant 
heat input of —117 BTU/(sq. ft.) (hr.). 

According to (29), if ¢ can be taken so small that all terms above the first are neg- 
lected, recrystallization cannot take place if k,¢’(0) > k.g(0). In fact, this condition 
would represent decrystallization and the curve would be in the negative x half plane. 


BOOK REVIEWS 


Electromagnetic theory. Proceedings of Symposia in Applied Mathematics, Volume II. 
American Mathematical Society, New York, 1950. iii + 91 pp. $3.00. 


This volume contains the proceedings of the Second Symposium in Applied Mathematics, held at 
the Massachusetts Institute of Technology in July 1948. Of the seventeen papers included, ten are pre- 
sented in the volume by abstract. 

The subject matter covers a wide range under the general topic of Electromagnetic Theory including 
related quantum theory and statistics. The papers in full are: H. Feshbach, The new quantum electrody- 
namics. J. L. Synge, Electromagnetism without metric. W. H. Watson, Discontinuity in electromagnetism. 
L. Infeld, The factorization method and its application to differential equations in theoretical physics. 
R. J. Duffin, Nonlinear electrical networks. C. L. Pekeris, Ray theory vs. normal mode theory in wave 
propagation problems. A. E. Heins, Systems of Wiener-Hopf integral equations and their application to 


some boundary value problems in electromagnetic theory. 


W. D. Hayes 
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Deformation and flow. An elementary introduction to theoretical rheology. By Markus 
teiner. H. K. Lewis & Co., Ltd., London, 1949. xix + 342 pp. 32s. 6d. (Sold in the 
U.S.A. by Interscience Publishers, New York, at $6.50). 


This treatment of the theory of Rheology is primarily concerned with the fundamental behavior of 
the various rheological media and not with the general formulations of the typical problems of the mechan- 
ics of continua. Thus, the kinematic considerations have been limited to deformations which are essen- 
tially pure shear but the range of media discussed is very broad. There are chapters dealing with viscous 
fluids, plastic media, the generalized Newtonian fluid, and complex media such as concrete; and there are 
chapters devoted to discussions of work hardening, creep, strength and failure, and the behavior near 
solid boundaries. The book contains a wealth of information for those interested in the behavior of various 
media under comparatively simple loading or simple flow conditions. 


G. F. CARRIER 


Experimental designs. By William G. Cochran and Gertrude M. Cox. John Wiley & Sons, 
Inc., New York and Chapman & Hall, Ltd., London. ix + 454 pp. $5.75. 


In this book the authors present a large number of different designs covering a great variety of experi- 
nental situations. The application of these designs and their analysis is illustrated by well chosen exam- 
ples. The book being intended as an elementary text and as a guide for research workers with little or no 
mathematical training, the authors do not fully discuss the basic assumptions underlying the analysis and 
the principles of statistical inference that lead from the assumptions to the analysis. However, the book 
distinguishes itself from many other elementary books in Statistics by giving a correct, although incom- 
plete discussion of the necessary fundamental statistical principles. 

The greatest value of the book in the reviewer’s opinion is the extensive collection of experimental 
designs presented. A disadvantage for the reader with more mathematical background is the absence of 
closed expressions for the computation of test functions. Such closed expressions, which would inform the 
mathematically trained reader at a glance, are replaced by numerical examples, probably because formu- 
lae involving multiple summation signs and other formidable mathematical symbols would deter the 
average reader. The reviewer believes, however, that the book will be very profitable reading also for the 
student of Mathematical Statistics, particularly when used as complementary reading to a treatise pre- 
senting the mathematical theory of the design of experiments and the analysis of variance. 


Henry B. MANN 


An introduction to the LaPlace transformation. By J. C. Jaeger. John Wiley & Sons, Inc., 
New York and Methuen & Co., Ltd., London. viii + 132 pp. $1.50. 


This member of the Methuen series of ““Monographs on physical subjects” is concerned with the use 
of the LaPlace transform in solving differential equations with constant coefficients. The first chapter 
introduces the theorems needed; the second deals with the analysis of problems associated with electric 
circuit analysis; the third introduces a few more theorems and the final chapter applies the operational 
technique to a partial differential equation. The concise presentation seems well suited to those who wish 
to become familiar with an operational tool without wading through a careful analysis of the pertinent 
mathematics, 


G. F. CARRIER 
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Elements of aerodynamics of supersonic flows. By Antonio Ferri. The Macmillan Company, 
New York, 1949. x + 434 pp. $10.00. 


This is a text designed for the practicing supersonic aerodynamicist. By limiting his field of coverage 
to purely supersonic flow, excepting incidentally subsonic flow where encountered behind a shock, the 
author has produced a treatment which may be termed reasonably comprehensive. The mathematical 
requirements for the reader are not too high—he needs a thorough knowledge of basic analysis through 
differential equations, but partial differential equations theory is not prerequisite and complex variables 
do not appear in the book. Although linearized theory is presented, the principal emphasis is on the non- 
linear, as befits a realistic approach. Test results are included wherever feasible for comparison and illus- 
tration and the entire presentation is made very graphic by a satisfying number of figures. The index is 
completely inadequate. 

This is definitely not a text for workers in fundamental research. The main reason for this is that the 
main fundamental problems of high speed flow lie outside the field of coverage of the book; they involve 
subsonic flow regions, viscous effects, and stability questions, and much of the mathematics involved is 
more advanced. Also, in designing his book for the engineer, the author has emphasized the specific at the 
expense of the fundamental, omitting such basic concepts in supersonic flow as those of the Euler mo- 
mentum integral and the linearized and hypersonic similitudes. 


W. D. Hayes 


Methoden der praktischen analysis. By F. A. Willers. 2nd ed. Walter de Gruyter & Co., 
Berlin, 1950. 410 pp. DM 24.-. 


The first. edition of this well-known text appeared in 1928. In preparing the present revised edition, 
the author has included many new methods developed since the publication of the first edition. To keep 
the size of the volume approximately the same, some special graphical methods discussed were omitted. 
While there are frequent minor revisions throughout the book, the major additions are sections on the 
approximate integration of integral equations of the second kind, evaluation of improper integrals, 
trigonometric interpolation, and the solution of systems of linear equations by iteration. The last chapter 
on approximate integration of differential equations has been completely rewritten. According to a state- 
ment in the preface, the manuscript of this revised edition was completed in 1943. This explains the fact 
that useful techniques developed during the war in the United States and Great Britain have not been 


included. 
W. PRAGER 


Tables of the generalized exponential-integral functions. By The Staff of the Computation 
Laboratory. Harvard University Press, Cambridge, Massachusetts, 1949. xxv + 416 
pp. $8.00. 


The twenty-five introductory pages of this volume are concerned with the definitions of the functions 
to be tabulated, a discussion of the computation methods, the interpolation technique and a brief applica- 
tion discussion with bibliography. Three functions are computed: They are {7F,(a, s) ds where the three 
integrands used are u~!(1 — e~), ute sin u, and u“1(1 — e™ cos u), with u(a, s) = (a? + s*)!/2, 

These functions are tabulated in the range 0 < a < .049 in steps Aa = .001 with O < xz < .049 in 
steps of .001;in0 < a < .098 and 0 < x < .098 with increments in x and a of .002; inO0 < a < .245 
and 0 < xz < .245 with increments of .005; for 0 < a < .49,0 < x < .49 with increments .01; in 
0 <a,z < .98 with increments .02;in 0 < a,z < 2.45 with increments .05; for 0 < a, x < .49 with 
increments .1; and for 0 < a, x < .98 with increments .2. 

G. F. Carrier 
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Tables of the function sin ¢/¢ and of its first eleven derivatives. By The Staff of the Computa- 
tion Laboratory. Harvard University Press, Cambridge, Massachusetts, 1949. sviii + 
241 pp. $8.00. 


The function sin ¢/@ and its first eleven derivatives are tabulated for values of ¢ in the range 


0 < @ < 20x. The increments in ¢ are taken to be .5 degrees for computation convenience. Remarks 


about computational procedure and a brief discussion of the applications in which these functions are 
useful make up the introductory pages. 
G. F. CARRIER 


Tables of the Bessel functions of the first kind of orders sixty-four through seventy-eight. By 
The Staff of the Computation Laboratory. Harvard University Press, Cambridge, 
Massachusetts, 1949. 566 pp. $8.00. 


The Bessel Functions J6s(x), Jes(x), «** Jzs(x) are tabulated to ten decimal places. The argument 
varies in steps of .0O1 forO < 2x < 100 
G. F. Carrier 


The variational principles of mechanics. By Cornelius Lanczos. University of Toronto 
Press, Toronto, 1949. xxv + 307 pp. $5.75. 


This is an excellent book covering one of the most brilliant achievements of scientific theory, and the 
text will certainly be welcomed by teachers and students of analytical mechanics. 

The author leads the reader through the historical development, starting from the basic concepts of 
mechanics and discussing a number of specific examples to familiarize the student with the general prin- 
ciples. Short, clear summaries present the fundamental results obtained in each section and make it 
possible to follow step by step the general line of thought. The chapter on Liouville’s theorem may be 
taken as a typical illustration of the author’s method, with a general discussion of phase-space, then a few 
simple examples, and finally the general theorem about integral-invariants. 

One may only regret that the author did not include Ehrenfest’s adiabatic invariants, which play 
such an important role in quantum theories. The analytical methods used in the old quantum mechanics 
are very completely and clearly discussed in the last chapter, which will be found especially valuable by 
many students in theoretical physics. 

This is a brilliant presentation of a great chapter of science, and the text will be of great help to 
students who have studied mechanics in Synge and Griffith’s excellent treatise where these problems were 


barely sketched. 


L. BRILLOUIN 


The meaning of relativity. By Albert Einstein. Princeton University Press, Princeton, New 
Jersey, 1950. 150 pp. $2.50. 


The first part of the book is a general explanation of the physical background of Relativity. It begins 
with a discussion of the fundamental notions of space and time in classical physics and the role of rela- 
tivity in classical mechanics, where all physical laws are assumed to be invariant for Cartesian systems of 
space coordinates, while the time variable is not supposed to be affected by these Cartesian changes of 


coordinates: at this stage, time and distance are absolute. 
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This set of definitions fails when applied to the Maxwell-Lorentz equations of electromagnetism. The 
discussion of this problem led Einstein to the discovery of the theory of Special Relativity, and the trend 
of thought is very clearly explained in the second chapter of the book, where the fusion of space and time 
in a comprehensive system of 4 variables is introduced. 

The famous relation E = me? is obtained, but there is just a short footnote to tell of the experimental 
proofs of this most important prediction. 

General relativity was developed in order to include gravitation in this general description. The point 
of departure is the identity of inert mass and gravitational mass. This leads directly to the principle of 
equivalence, which states that gravitation plays a similar role to acceleration, and explains the need for 

tiemannian geometry. The connection between the General Relativity and Mach’s ideas makes for a 
very interesting discussion of the fundamental assumptions and of the new results to be expected. The 
chapter concludes in favor of a space-bounded universe, as representing the most satisfactory type of 
general equation, 

In the appendix to the second edition, this cosmologic problem is discussed, in connection with the 
experimental results about the red shift of spectral lines from distant nebulae and the expanding universe. 
The author does not believe in the addition of a “cosmologic term” in the equations of gravity, and pre- 
sents the theory of Friedman as the best solution, with a four-dimensional universe that is isotropic with 
respect to three dimensions (space). This yields a value of 10° years for the age of the universe, a very 
small value indeed and difficult to reconcile with other astronomical data. 

The second appendix is a new addition to the present third edition, and contains the new generalized 
theory of gravitation, recently developed by A. Einstein. This is an attempt to include Electromagnetism 
in the general geometry of the universe, instead of leaving Maxwell’s equations as a separate group of 
equations, unconnected with the geometry. The previous general Relativity was a frame in which Electro- 
magnetism was included as a separate picture. The aim is now to rebuild frame and picture together as ¢ 
single unit. For this problem, our general knowledge of physics does not permit an unequivocal choice as 
did the principle of equivalence for gravitation. The only clue is that something similar to Maxwell’s 
equation must be contained within the total field. 

Many such attempts have been presented in the last decades, and they can be classified into two main 
groups. The electromagnetic field representing a skew-symmetrical tensor, one may try to add such skew- 
symmetrical terms to either 

A—the g,, fundamental metric tensor 


or B—the T‘,, coefficients of parallel displacement 


uy 


(the T',, do not constitute a tensor) 


Both quantities were usually defined as symmetrical in wv, and can be completed with skew-symmetrical 
terms. Einstein uses both generalizations simultaneously and builds up the corresponding Riemann Ten- 
sor for curvature. He then looks for field equations, and in order to be certain of the compatibility of these 
equations he derives them from a variational principle according to the well-known Hamilton procedure. 
After some simplification, the result is a system of equations very similar to Maxwell’s equations (with no 
free magnetic charges). 

The author does not give any discussion of the practical value of these new equations and does not 
tell the new facts that they may include. One must remember that the theoretical problem was a formid- 
able one, and that A. Einstein has been working on it for two decades, trying all possible lines of attack 
and attempting to use all sorts of new geometrical structures. He finally makes a choice, and gives us a set 
of formulas that seem to him to represent the best logical solution, leaving to others the trouble to discuss 
the practical and physical meaning of the new theory. Many questions remain still obscure, and first of 
all: how is it possible to obtain the equations of motion for different particles with different masses, 
charges and magnetic moments? One single geodesic can not do, since all trajectories are different from 
each other. This problem is not discussed in Einstein’s paper, and we hope he shall give us his ideas about 
the type of logical solution to be expected. 

The reader will refer with great interest to the paper published by Einstein in the “Scientific Amer- 
ican’’, p. 13, vol. 182, no. 4, April 1950, where the philosophical background of the theory is discussed with 
great care and remarkable clarity. Einstein indicates in this paper a very curious feature of his theory: the 
initial conditions cannot be freely chosen. He points out that this fact does not seem to lead to serious 
trouble, provided one may prove that the system is extensive enough to cover all the possible physical 


problems, and this is still an open question. 
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Let us admire the genius of a great thinker and an exceptional scientist, who fulfills now, in his ripe 
age, an audacious plan of research that he started as a young man. A wonderful achievement and a 
rare one, 


L. BRILLOUIN 


Physical principles of oil production. By Morris Muskat. First Edition. McGraw-Hill 
Book Company, Inc., New York, Toronto, London, 1949. xv + 922 pp. $15.00. 


As stated in the preface this work is neither a textbook nor a reorganized version of well established 


theories for pedagogical purposes, but an exposition of material which is in an active state of flux, with the 
purpose of formulating and correlating the physical principles and facts underlying the mechanisms of oil 


production. 

The physics of oil production has been a rapidly expanding branch of applied physics in the last 
fifteen years and one cannot help but admire the art and ingenuity with which the scientific approach has 
been introduced in a field where rule of thumb used to prevail and where the object of study, the oil reser- 
voir, offers an infinite variety of individual aspects which are susceptible to only limited observation. 

Chapters 1 to 3 describe and discuss the general physical properties of petroleum fluids and oil bearing 
rocks with particular reference to the thermodynamic behavior of multiphase hydrocarbon systems and 
their interaction with porous rock. The following three chapters deal with the motion of single phase fluids 
through porous rock. Incompressible fluid flows form the most important class for practical applications 
and such flows in homogeneous and heterogeneous media are treated by solving the Laplace equation. 
Many types of flow fields are described corresponding to various physical conditions of stratification and 
well penetration. The effect of compressibility is also introduced and the corresponding diffusion-type 
equation is applied to a few idealized cases. The more difficult problems of mult phase fluid flow phenom- 
ena are treated in Chapters 7 to 11. The physical concepts and characteristics of this type of flow are 
reviewed and applied to the solution of several simple steady-state flow equations. This is supplemented 
by examples of actual performance records of fields. The last three chapters are devoted to problems of 
secondary recovery, condensate reservoirs, well spacing and the estimation of recoverable reserves. The 
mathematical and analogue computer methods developed for the solution of these problems will be of 
interest also to those who are not directly concerned with oil production. 

In view of the authoritative nature of this work it is unfortunate that serious criticism may be ex- 
pressed concerning the presentation of the general laws of fluid flow through porous media. The treatment 
lacks the clarity and generality which one expects to find in a treatise of this scope. In the opinion of the 
reviewer the book would have gained appreciably if it had followed a more systematic presentation as 
developed, e.g., by M. K. Hubbert (cf. Theory of Ground Water Motion, The Journal of Geology Vol. 48, 
No. 8, November-December 1940). 

A great deal of interest has been attracted lately to the question of the rapid depletion of the world 
oil reserves. The challenge to the applied mathematician and physicist is well illustrated by this most 
thoroughly documented treatment of the field. 

M. A. Bior 
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